Changeset 4648 for palm


Ignore:
Timestamp:
Aug 25, 2020 7:52:08 AM (2 months ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

Location:
palm/trunk/SOURCE
Files:
11 edited

Legend:

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

    r4548 r4648  
    11!> @file init_3d_model.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    2120! ------------------
    22 ! 
    23 ! 
     21!
     22!
    2423! Former revisions:
    2524! -----------------
    2625! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4548 2020-05-28 19:36:45Z suehring
    2729! Bugfix, move call for lsf_forcing_surf after lsf_init is called
    28 ! 
     30!
    2931! 4514 2020-04-30 16:29:59Z suehring
    30 ! Add possibility to initialize surface sensible and latent heat fluxes via
    31 ! a static driver.
    32 !
     32! Add possibility to initialize surface sensible and latent heat fluxes via a static driver.
     33!
    3334! 4493 2020-04-10 09:49:43Z pavelkrc
    34 ! Overwrite u_init, v_init, pt_init, q_init and s_init with hom for all
    35 ! cyclic_fill-cases, not only for turbulent_inflow = .TRUE.
    36 ! 
     35! Overwrite u_init, v_init, pt_init, q_init and s_init with hom for all cyclic_fill-cases, not only
     36! for turbulent_inflow = .TRUE.
     37!
    3738! 4360 2020-01-07 11:25:50Z suehring
    38 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    39 ! topography information used in wall_flags_static_0
    40 ! 
     39! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     40! information used in wall_flags_static_0
     41!
    4142! 4329 2019-12-10 15:46:36Z motisi
    4243! Renamed wall_flags_0 to wall_flags_static_0
    43 ! 
     44!
    4445! 4286 2019-10-30 16:01:14Z resler
    4546! implement new palm_date_time_mod
    46 ! 
     47!
    4748! 4223 2019-09-10 09:20:47Z gronemeier
    48 ! Deallocate temporary string array since it may be re-used to read different
    49 ! input data in other modules
    50 ! 
     49! Deallocate temporary string array since it may be re-used to read different input data in other
     50! modules
     51!
    5152! 4186 2019-08-23 16:06:14Z suehring
    52 ! Design change, use variables defined in netcdf_data_input_mod to read netcd
    53 ! variables rather than define local ones.
    54 ! 
     53! Design change, use variables defined in netcdf_data_input_mod to read netcd variables rather than
     54! define local ones.
     55!
    5556! 4185 2019-08-23 13:49:38Z oliver.maas
    5657! For initializing_actions = ' cyclic_fill':
    57 ! Overwrite u_init, v_init, pt_init, q_init and s_init with the
    58 ! (temporally) and horizontally averaged vertical profiles from the end
    59 ! of the prerun, because these profiles shall be used as the basic state
    60 ! for the rayleigh damping and the pt_damping.
    61 !
     58! Overwrite u_init, v_init, pt_init, q_init and s_init with the (temporally) and horizontally
     59! averaged vertical profiles from the end of the prerun, because these profiles shall be used as the
     60! basic state for the rayleigh damping and the pt_damping.
     61!
    6262! 4182 2019-08-22 15:20:23Z scharf
    6363! Corrected "Former revisions" section
    64 ! 
     64!
    6565! 4168 2019-08-16 13:50:17Z suehring
    6666! Replace function get_topography_top_index by topo_top_ind
    67 ! 
     67!
    6868! 4151 2019-08-09 08:24:30Z suehring
    6969! Add netcdf directive around input calls (fix for last commit)
    70 ! 
     70!
    7171! 4150 2019-08-08 20:00:47Z suehring
    72 ! Input of additional surface variables independent on land- or urban-surface
    73 ! model
    74 !
     72! Input of additional surface variables independent on land- or urban-surface model
     73!
    7574! 4131 2019-08-02 11:06:18Z monakurppa
    7675! Allocate sums and sums_l to allow profile output for salsa variables.
    77 ! 
     76!
    7877! 4130 2019-08-01 13:04:13Z suehring
    79 ! Effectively reduce 3D initialization to 1D initial profiles. This is because
    80 ! 3D initialization produces structures in the w-component that are correlated
    81 ! with the processor grid for some unknown reason 
    82 ! 
     78! Effectively reduce 3D initialization to 1D initial profiles. This is because 3D initialization
     79! produces structures in the w-component that are correlated with the processor grid for some
     80! unknown reason
     81!
    8382! 4090 2019-07-11 15:06:47Z Giersch
    8483! Unused variables removed
    85 ! 
     84!
    8685! 4088 2019-07-11 13:57:56Z Giersch
    8786! Pressure and density profile calculations revised using basic functions
    88 ! 
     87!
    8988! 4048 2019-06-21 21:00:21Z knoop
    9089! Further modularization of particle code components
    91 ! 
     90!
    9291! 4017 2019-06-06 12:16:46Z schwenkel
    93 ! Convert most location messages to debug messages to reduce output in
    94 ! job logfile to a minimum
    95 !
    96 !
     92! Convert most location messages to debug messages to reduce output in job logfile to a minimum
     93!
    9794! unused variable removed
    98 ! 
     95!
    9996! 3937 2019-04-29 15:09:07Z suehring
    100 ! Move initialization of synthetic turbulence generator behind initialization
    101 ! of offline nesting. Remove call for stg_adjust, as this is now already done
    102 ! in stg_init.
    103 !
     97! Move initialization of synthetic turbulence generator behind initialization of offline nesting.
     98! Remove call for stg_adjust, as this is now already done in stg_init.
     99!
    104100! 3900 2019-04-16 15:17:43Z suehring
    105101! Fix problem with LOD = 2 initialization
    106 ! 
     102!
    107103! 3885 2019-04-11 11:29:34Z kanani
    108 ! Changes related to global restructuring of location messages and introduction
    109 ! of additional debug messages
    110 ! 
     104! Changes related to global restructuring of location messages and introduction of additional debug
     105! messages
     106!
    111107! 3849 2019-04-01 16:35:16Z knoop
    112108! Move initialization of rmask before initializing user_init_arrays
    113 ! 
     109!
    114110! 3711 2019-01-31 13:44:26Z knoop
    115111! Introduced module_interface_init_checks for post-init checks in modules
    116 ! 
     112!
    117113! 3700 2019-01-26 17:03:42Z knoop
    118114! Some interface calls moved to module_interface + cleanup
    119 ! 
     115!
    120116! 3648 2019-01-02 16:35:46Z suehring
    121117! Rename subroutines for surface-data output
     
    133129!> or
    134130!> c) read values of a previous run
    135 !------------------------------------------------------------------------------!
     131!--------------------------------------------------------------------------------------------------!
    136132 SUBROUTINE init_3d_model
    137133
     
    141137    USE arrays_3d
    142138
    143     USE basic_constants_and_equations_mod,                                     &
    144         ONLY:  c_p, g, l_v, pi, exner_function, exner_function_invers,         &
    145                ideal_gas_law_rho, ideal_gas_law_rho_pt, barometric_formula
    146 
    147     USE bulk_cloud_model_mod,                                                  &
     139    USE basic_constants_and_equations_mod,                                                         &
     140        ONLY:  barometric_formula, c_p, exner_function, exner_function_invers, g,                  &
     141               ideal_gas_law_rho, ideal_gas_law_rho_pt, l_v, pi
     142
     143    USE bulk_cloud_model_mod,                                                                      &
    148144        ONLY:  bulk_cloud_model
    149145
    150     USE chem_modules,                                                          &
     146    USE chem_modules,                                                                              &
    151147        ONLY:  max_pr_cs ! ToDo: this dependency needs to be removed cause it is ugly #new_dom
    152148
    153149    USE control_parameters
    154150
    155     USE grid_variables,                                                        &
     151    USE grid_variables,                                                                            &
    156152        ONLY:  dx, dy, ddx2_mg, ddy2_mg
    157153
     
    159155
    160156    USE kinds
    161  
    162     USE lsf_nudging_mod,                                                       &
     157
     158    USE lsf_nudging_mod,                                                                           &
    163159        ONLY:  ls_forcing_surf
    164160
    165     USE model_1d_mod,                                                          &
     161    USE model_1d_mod,                                                                              &
    166162        ONLY:  init_1d_model, l1d, u1d, v1d
    167163
    168     USE module_interface,                                                      &
    169         ONLY:  module_interface_init_arrays,                                   &
    170                module_interface_init,                                          &
     164    USE module_interface,                                                                          &
     165        ONLY:  module_interface_init_arrays,                                                       &
     166               module_interface_init,                                                              &
    171167               module_interface_init_checks
    172168
    173     USE multi_agent_system_mod,                                                &
     169    USE multi_agent_system_mod,                                                                    &
    174170        ONLY:  agents_active, mas_init
    175171
    176     USE netcdf_interface,                                                      &
     172    USE netcdf_interface,                                                                          &
    177173        ONLY:  dots_max
    178174
    179     USE netcdf_data_input_mod,                                                 &
    180         ONLY:  char_fill,                                                      &
    181                check_existence,                                                &
    182                close_input_file,                                               &
    183                get_attribute,                                                  &
    184                get_variable,                                                   &
    185                init_3d,                                                        &
    186                input_pids_static,                                              &
    187                inquire_num_variables,                                          &
    188                inquire_variable_names,                                         &
    189                input_file_static,                                              &
    190                netcdf_data_input_init_3d,                                      &
    191                num_var_pids,                                                   &
    192                open_read_file,                                                 &
    193                pids_id,                                                        &
    194                real_2d,                                                        &
     175    USE netcdf_data_input_mod,                                                                     &
     176        ONLY:  char_fill,                                                                          &
     177               check_existence,                                                                    &
     178               close_input_file,                                                                   &
     179               get_attribute,                                                                      &
     180               get_variable,                                                                       &
     181               init_3d,                                                                            &
     182               input_pids_static,                                                                  &
     183               inquire_num_variables,                                                              &
     184               inquire_variable_names,                                                             &
     185               input_file_static,                                                                  &
     186               netcdf_data_input_init_3d,                                                          &
     187               num_var_pids,                                                                       &
     188               open_read_file,                                                                     &
     189               pids_id,                                                                            &
     190               real_2d,                                                                            &
    195191               vars_pids
    196                
    197     USE nesting_offl_mod,                                                      &
     192
     193    USE nesting_offl_mod,                                                                          &
    198194        ONLY:  nesting_offl_init
    199195
    200     USE palm_date_time_mod,                                                    &
     196    USE palm_date_time_mod,                                                                        &
    201197        ONLY:  set_reference_date_time
    202198
     
    204200
    205201#if defined( __parallel )
    206     USE pmc_interface,                                                         &
     202    USE pmc_interface,                                                                             &
    207203        ONLY:  nested_run
    208204#endif
    209205
    210     USE random_function_mod 
    211 
    212     USE random_generator_parallel,                                             &
     206    USE random_function_mod
     207
     208    USE random_generator_parallel,                                                                 &
    213209        ONLY:  init_parallel_random_generator
    214210
    215     USE read_restart_data_mod,                                                 &
    216         ONLY:  rrd_read_parts_of_global, rrd_local
    217 
    218     USE statistics,                                                            &
    219         ONLY:  hom, hom_sum, mean_surface_level_height, pr_palm, rmask,        &
    220                statistic_regions, sums, sums_divnew_l, sums_divold_l, sums_l,  &
    221                sums_l_l, sums_wsts_bc_l, ts_value,                             &
     211    USE read_restart_data_mod,                                                                     &
     212        ONLY:  rrd_local, rrd_read_parts_of_global
     213
     214    USE statistics,                                                                                &
     215        ONLY:  hom, hom_sum, mean_surface_level_height, pr_palm, rmask, statistic_regions, sums,   &
     216               sums_divnew_l, sums_divold_l, sums_l, sums_l_l, sums_wsts_bc_l, ts_value,           &
    222217               weight_pres, weight_substep
    223218
    224     USE synthetic_turbulence_generator_mod,                                    &
     219    USE synthetic_turbulence_generator_mod,                                                        &
    225220        ONLY:  stg_init, use_syn_turb_gen
    226221
    227     USE surface_layer_fluxes_mod,                                              &
     222    USE surface_layer_fluxes_mod,                                                                  &
    228223        ONLY:  init_surface_layer_fluxes
    229224
    230     USE surface_mod,                                                           &
    231         ONLY :  init_single_surface_properties,                                &
    232                 init_surface_arrays,                                           &
    233                 init_surfaces,                                                 &
    234                 surf_def_h,                                                    &
    235                 surf_def_v,                                                    &
    236                 surf_lsm_h,                                                    &
     225    USE surface_mod,                                                                               &
     226        ONLY :  init_single_surface_properties,                                                    &
     227                init_surface_arrays,                                                               &
     228                init_surfaces,                                                                     &
     229                surf_def_h,                                                                        &
     230                surf_def_v,                                                                        &
     231                surf_lsm_h,                                                                        &
    237232                surf_usm_h
    238233
    239234#if defined( _OPENACC )
    240     USE surface_mod,                                                           &
     235    USE surface_mod,                                                                               &
    241236        ONLY :  bc_h
    242237#endif
    243238
    244     USE surface_data_output_mod,                                               &
     239    USE surface_data_output_mod,                                                                   &
    245240        ONLY:  surface_data_output_init
    246241
    247242    USE transpose_indices
    248243
     244
    249245    IMPLICIT NONE
    250    
     246
    251247    INTEGER(iwp) ::  i                    !< grid index in x direction
    252248    INTEGER(iwp) ::  ind_array(1)         !< dummy used to determine start index for external pressure forcing
     
    254250    INTEGER(iwp) ::  k                    !< grid index in z direction
    255251    INTEGER(iwp) ::  k_surf               !< surface level index
    256     INTEGER(iwp) ::  l                    !< running index over surface orientation   
    257     INTEGER(iwp) ::  m                    !< index of surface element in surface data type   
     252    INTEGER(iwp) ::  l                    !< running index over surface orientation
     253    INTEGER(iwp) ::  m                    !< index of surface element in surface data type
    258254    INTEGER(iwp) ::  nz_u_shift           !< topography-top index on u-grid, used to vertically shift initial profiles
    259255    INTEGER(iwp) ::  nz_v_shift           !< topography-top index on v-grid, used to vertically shift initial profiles
     
    267263    INTEGER(iwp) ::  sr                   !< index of statistic region
    268264
    269     INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  ngp_2dh_l  !< toal number of horizontal grid points in statistical region on subdomain
    270 
    271     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l    !< number of horizontal non-wall bounded grid points on subdomain
    272     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_s_inner_l  !< number of horizontal non-topography grid points on subdomain
    273 
    274 
    275    
    276     REAL(wp), DIMENSION(:), ALLOCATABLE ::  init_l        !< dummy array used for averaging 3D data to obtain inital profiles
    277     REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_hydrostatic !< hydrostatic pressure
     265    INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  ngp_2dh_l  !< toal number of horizontal grid points in statistical region on
     266                                                             !< subdomain
     267
     268    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l    !< number of horizontal non-wall bounded grid points on
     269                                                                     !< subdomain
     270    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_s_inner_l  !< number of horizontal non-topography grid points on
     271                                                                     !< subdomain
    278272
    279273    REAL(wp) ::  dx_l !< grid spacing along x on different multigrid level
    280274    REAL(wp) ::  dy_l !< grid spacing along y on different multigrid level
    281275
     276    REAL(wp), DIMENSION(:), ALLOCATABLE ::  init_l                       !< dummy array used for averaging 3D data to obtain
     277                                                                         !< inital profiles
     278    REAL(wp), DIMENSION(:), ALLOCATABLE ::  mean_surface_level_height_l  !< mean surface level height on subdomain
     279    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_l               !< total number of non-topography grid points on subdomain
     280    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_tmp             !< total number of non-topography grid points
     281    REAL(wp), DIMENSION(:), ALLOCATABLE ::  p_hydrostatic                !< hydrostatic pressure
     282
    282283    REAL(wp), DIMENSION(1:3) ::  volume_flow_area_l     !< area of lateral and top model domain surface on local subdomain
    283284    REAL(wp), DIMENSION(1:3) ::  volume_flow_initial_l  !< initial volume flow into model domain
    284285
    285     REAL(wp), DIMENSION(:), ALLOCATABLE ::  mean_surface_level_height_l !< mean surface level height on subdomain
    286     REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_l    !< total number of non-topography grid points on subdomain
    287     REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_tmp  !< total number of non-topography grid points
    288 
    289286    TYPE(real_2d) ::  tmp_2d !< temporary variable to input additional surface-data from static file
    290    
     287
    291288    CALL location_message( 'model initialization', 'start' )
    292289!
     
    297294!
    298295!-- Allocate arrays
    299     ALLOCATE( mean_surface_level_height(0:statistic_regions),                  &
    300               mean_surface_level_height_l(0:statistic_regions),                &
    301               ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions),    &
    302               ngp_3d(0:statistic_regions),                                     &
    303               ngp_3d_inner(0:statistic_regions),                               &
    304               ngp_3d_inner_l(0:statistic_regions),                             &
    305               ngp_3d_inner_tmp(0:statistic_regions),                           &
    306               sums_divnew_l(0:statistic_regions),                              &
     296    ALLOCATE( mean_surface_level_height(0:statistic_regions),                                      &
     297              mean_surface_level_height_l(0:statistic_regions),                                    &
     298              ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions),                        &
     299              ngp_3d(0:statistic_regions),                                                         &
     300              ngp_3d_inner(0:statistic_regions),                                                   &
     301              ngp_3d_inner_l(0:statistic_regions),                                                 &
     302              ngp_3d_inner_tmp(0:statistic_regions),                                               &
     303              sums_divnew_l(0:statistic_regions),                                                  &
    307304              sums_divold_l(0:statistic_regions) )
    308305    ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt) )
    309     ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                    &
    310               ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),                  &
    311               ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),                  &
    312               ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),                &
    313               rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),                  &
    314               sums(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa),      &
     306    ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                                        &
     307              ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),                                      &
     308              ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),                                      &
     309              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),                                    &
     310              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),                                      &
     311              sums(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa),                          &
    315312              sums_l(nzb:nzt+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0:threads_per_task-1),   &
    316               sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1),    &
     313              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1),                        &
    317314              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions) )
    318315    ALLOCATE( ts_value(dots_max,0:statistic_regions) )
    319316    ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) )
    320317
    321     ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr),                                    &
    322               p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
     318    ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr),                                                        &
     319              p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                    &
    323320              tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    324321
    325     ALLOCATE( pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
    326               pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
    327               u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
    328               u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
    329               u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
    330               v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
    331               v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
    332               v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
    333               w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
    334               w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
     322    ALLOCATE( pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                 &
     323              pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                 &
     324              u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
     325              u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
     326              u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
     327              v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
     328              v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
     329              v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
     330              w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
     331              w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                  &
    335332              w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    336333    IF (  .NOT.  neutral )  THEN
     
    339336!
    340337!-- Pre-set masks for regional statistics. Default is the total model domain.
    341 !-- Ghost points are excluded because counting values at the ghost boundaries
    342 !-- would bias the statistics
     338!-- Ghost points are excluded because counting values at the ghost boundaries would bias the
     339!-- statistics.
    343340    rmask = 1.0_wp
    344341    rmask(:,nxlg:nxl-1,:) = 0.0_wp;  rmask(:,nxr+1:nxrg,:) = 0.0_wp
    345342    rmask(nysg:nys-1,:,:) = 0.0_wp;  rmask(nyn+1:nyng,:,:) = 0.0_wp
    346343!
    347 !-- Following array is required for perturbation pressure within the iterative
    348 !-- pressure solvers. For the multistep schemes (Runge-Kutta), array p holds
    349 !-- the weighted average of the substeps and cannot be used in the Poisson
    350 !-- solver.
     344!-- Following array is required for perturbation pressure within the iterative pressure solvers. For
     345!-- the multistep schemes (Runge-Kutta), array p holds the weighted average of the substeps and
     346!-- cannot be used in the Poisson solver.
    351347    IF ( psolver == 'sor' )  THEN
    352348       ALLOCATE( p_loc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     
    367363!
    368364!--    3D-humidity
    369        ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
    370                  q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
    371                  q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
    372                  vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 
    373     ENDIF   
    374    
     365       ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                               &
     366                 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                               &
     367                 q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                               &
     368                 vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     369    ENDIF
     370
    375371    IF ( passive_scalar )  THEN
    376372
    377373!
    378374!--    3D scalar arrays
    379        ALLOCATE( s_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
    380                  s_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
     375       ALLOCATE( s_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                               &
     376                 s_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                               &
    381377                 s_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    382378
     
    384380
    385381!
    386 !-- Allocate and set 1d-profiles for Stokes drift velocity. It may be set to
    387 !-- non-zero values later in ocean_init
    388     ALLOCATE( u_stokes_zu(nzb:nzt+1), u_stokes_zw(nzb:nzt+1),                  &
     382!-- Allocate and set 1d-profiles for Stokes drift velocity. It may be set to non-zero values later
     383!-- in ocean_init.
     384    ALLOCATE( u_stokes_zu(nzb:nzt+1), u_stokes_zw(nzb:nzt+1),                                      &
    389385              v_stokes_zu(nzb:nzt+1), v_stokes_zw(nzb:nzt+1) )
    390386    u_stokes_zu(:) = 0.0_wp
     
    401397    ALLOCATE( drho_air_zw(nzb:nzt+1) )
    402398!
    403 !-- Density profile calculation for anelastic and Boussinesq approximation
    404 !-- In case of a Boussinesq approximation, a constant density is calculated
    405 !-- mainly for output purposes. This density do not need to be considered
    406 !-- in the model's system of equations.
     399!-- Density profile calculation for anelastic and Boussinesq approximation.
     400!-- In case of a Boussinesq approximation, a constant density is calculated mainly for output
     401!-- purposes. This density does not need to be considered in the model's system of equations.
    407402    IF ( TRIM( approximation ) == 'anelastic' )  THEN
    408403       DO  k = nzb, nzt+1
    409           p_hydrostatic(k) = barometric_formula(zu(k), pt_surface *            &
    410                              exner_function(surface_pressure * 100.0_wp),      &
    411                              surface_pressure * 100.0_wp)
    412          
     404          p_hydrostatic(k) = barometric_formula(zu(k), pt_surface *                                &
     405                                                exner_function(surface_pressure * 100.0_wp),       &
     406                                                surface_pressure * 100.0_wp)
     407
    413408          rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), pt_init(k))
    414409       ENDDO
    415        
     410
    416411       DO  k = nzb, nzt
    417412          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
    418413       ENDDO
    419        
    420        rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
    421                             + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
    422                            
     414
     415       rho_air_zw(nzt+1)  = rho_air_zw(nzt) + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
     416
    423417    ELSE
    424418       DO  k = nzb, nzt+1
    425           p_hydrostatic(k) = barometric_formula(zu(nzb), pt_surface *          &
    426                              exner_function(surface_pressure * 100.0_wp),      &
    427                              surface_pressure * 100.0_wp)
     419          p_hydrostatic(k) = barometric_formula(zu(nzb), pt_surface *                              &
     420                                                exner_function(surface_pressure * 100.0_wp),       &
     421                                                surface_pressure * 100.0_wp)
    428422
    429423          rho_air(k) = ideal_gas_law_rho_pt(p_hydrostatic(k), pt_init(nzb))
    430424       ENDDO
    431        
     425
    432426       DO  k = nzb, nzt
    433427          rho_air_zw(k) = 0.5_wp * ( rho_air(k) + rho_air(k+1) )
    434428       ENDDO
    435        
    436        rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
    437                             + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
    438                            
    439     ENDIF
    440 !
    441 !-- compute the inverse density array in order to avoid expencive divisions
     429
     430       rho_air_zw(nzt+1)  = rho_air_zw(nzt) + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
     431
     432    ENDIF
     433!
     434!-- Compute the inverse density array in order to avoid expencive divisions
    442435    drho_air    = 1.0_wp / rho_air
    443436    drho_air_zw = 1.0_wp / rho_air_zw
     
    453446
    454447!
    455 !-- calculate flux conversion factors according to approximation and in-/output mode
     448!-- Calculate flux conversion factors according to approximation and in-/output mode
    456449    DO  k = nzb, nzt+1
    457450
     
    484477
    485478!
    486 !-- In case of multigrid method, compute grid lengths and grid factors for the
    487 !-- grid levels with respective density on each grid
     479!-- In case of multigrid method, compute grid lengths and grid factors for the grid levels with
     480!-- respective density on each grid.
    488481    IF ( psolver(1:9) == 'multigrid' )  THEN
    489482
     
    500493       dzu_mg(:,maximum_grid_level) = dzu
    501494       rho_air_mg(:,maximum_grid_level) = rho_air
    502 !       
    503 !--    Next line to ensure an equally spaced grid. 
     495!
     496!--    Next line to ensure an equally spaced grid.
    504497       dzu_mg(1,maximum_grid_level) = dzu(2)
    505        rho_air_mg(nzb,maximum_grid_level) = rho_air(nzb) +                     &
    506                                              (rho_air(nzb) - rho_air(nzb+1))
     498       rho_air_mg(nzb,maximum_grid_level) = rho_air(nzb) + (rho_air(nzb) - rho_air(nzb+1))
    507499
    508500       dzw_mg(:,maximum_grid_level) = dzw
     
    512504           dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1)
    513505           dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1)
    514            rho_air_mg(nzb,l)    = rho_air_mg(nzb,l+1) + (rho_air_mg(nzb,l+1) - rho_air_mg(nzb+1,l+1))
    515            rho_air_zw_mg(nzb,l) = rho_air_zw_mg(nzb,l+1) + (rho_air_zw_mg(nzb,l+1) - rho_air_zw_mg(nzb+1,l+1))
     506           rho_air_mg(nzb,l)    = rho_air_mg(nzb,l+1)    + ( rho_air_mg(nzb,l+1)    -              &
     507                                                             rho_air_mg(nzb+1,l+1)    )
     508           rho_air_zw_mg(nzb,l) = rho_air_zw_mg(nzb,l+1) + ( rho_air_zw_mg(nzb,l+1) -              &
     509                                                             rho_air_zw_mg(nzb+1,l+1) )
    516510           rho_air_mg(nzb+1,l)    = rho_air_mg(nzb+1,l+1)
    517511           rho_air_zw_mg(nzb+1,l) = rho_air_zw_mg(nzb+1,l+1)
     
    534528             f2_mg(k,l) = rho_air_zw_mg(k,l) / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
    535529             f3_mg(k,l) = rho_air_zw_mg(k-1,l) / ( dzu_mg(k,l)   * dzw_mg(k,l) )
    536              f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) &
     530             f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) )                                     &
    537531                          * rho_air_mg(k,l) + f2_mg(k,l) + f3_mg(k,l)
    538532          ENDDO
     
    552546
    553547!
    554 !-- Arrays to store velocity data from t-dt and the phase speeds which
    555 !-- are needed for radiation boundary conditions
     548!-- Arrays to store velocity data from t-dt and the phase speeds which are needed for radiation
     549!-- boundary conditions.
    556550    IF ( bc_radiation_l )  THEN
    557        ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2),                               &
    558                  v_m_l(nzb:nzt+1,nysg:nyng,0:1),                               &
     551       ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2),                                                   &
     552                 v_m_l(nzb:nzt+1,nysg:nyng,0:1),                                                   &
    559553                 w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
    560554    ENDIF
    561555    IF ( bc_radiation_r )  THEN
    562        ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
    563                  v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
     556       ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                                               &
     557                 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                                               &
    564558                 w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
    565559    ENDIF
    566560    IF ( bc_radiation_l  .OR.  bc_radiation_r )  THEN
    567        ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng),           &
    568                  c_w(nzb:nzt+1,nysg:nyng) )
     561       ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng), c_w(nzb:nzt+1,nysg:nyng) )
    569562    ENDIF
    570563    IF ( bc_radiation_s )  THEN
    571        ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg),                               &
    572                  v_m_s(nzb:nzt+1,1:2,nxlg:nxrg),                               &
     564       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg),                                                   &
     565                 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg),                                                   &
    573566                 w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
    574567    ENDIF
    575568    IF ( bc_radiation_n )  THEN
    576        ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
    577                  v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
     569       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                                               &
     570                 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                                               &
    578571                 w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
    579572    ENDIF
    580573    IF ( bc_radiation_s  .OR.  bc_radiation_n )  THEN
    581        ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg),           &
    582                  c_w(nzb:nzt+1,nxlg:nxrg) )
    583     ENDIF
    584     IF ( bc_radiation_l  .OR.  bc_radiation_r  .OR.  bc_radiation_s  .OR.      &
    585          bc_radiation_n )  THEN
    586        ALLOCATE( c_u_m_l(nzb:nzt+1), c_v_m_l(nzb:nzt+1), c_w_m_l(nzb:nzt+1) )                   
     574       ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg), c_w(nzb:nzt+1,nxlg:nxrg) )
     575    ENDIF
     576    IF ( bc_radiation_l  .OR.  bc_radiation_r  .OR.  bc_radiation_s  .OR.  bc_radiation_n )  THEN
     577       ALLOCATE( c_u_m_l(nzb:nzt+1), c_v_m_l(nzb:nzt+1), c_w_m_l(nzb:nzt+1) )
    587578       ALLOCATE( c_u_m(nzb:nzt+1), c_v_m(nzb:nzt+1), c_w_m(nzb:nzt+1) )
    588579    ENDIF
     
    603594       vpt  => vpt_1
    604595    ENDIF
    605    
     596
    606597    IF ( passive_scalar )  THEN
    607598       s => s_1;  s_p => s_2;  ts_m => s_3
    608     ENDIF   
     599    ENDIF
    609600
    610601!
     
    617608
    618609!
    619 !-- Allocate arrays containing the RK coefficient for calculation of
    620 !-- perturbation pressure and turbulent fluxes. At this point values are
    621 !-- set for pressure calculation during initialization (where no timestep
    622 !-- is done). Further below the values needed within the timestep scheme
    623 !-- will be set.
    624     ALLOCATE( weight_substep(1:intermediate_timestep_count_max),               &
     610!-- Allocate arrays containing the RK coefficient for calculation of perturbation pressure and
     611!-- turbulent fluxes. At this point values are set for pressure calculation during initialization
     612!-- (where no timestep is done). Further below the values needed within the timestep scheme will be
     613!-- set.
     614    ALLOCATE( weight_substep(1:intermediate_timestep_count_max),                                   &
    625615              weight_pres(1:intermediate_timestep_count_max) )
    626616    weight_substep = 1.0_wp
    627617    weight_pres    = 1.0_wp
    628618    intermediate_timestep_count = 0  ! needed when simulated_time = 0.0
    629        
     619
    630620    IF ( debug_output )  CALL debug_message( 'allocating arrays', 'end' )
    631621
     
    636626!
    637627!-- Initialize local summation arrays for routine flow_statistics.
    638 !-- This is necessary because they may not yet have been initialized when they
    639 !-- are called from flow_statistics (or - depending on the chosen model run -
    640 !-- are never initialized)
    641     sums_divnew_l      = 0.0_wp
    642     sums_divold_l      = 0.0_wp
    643     sums_l_l           = 0.0_wp
    644     sums_wsts_bc_l     = 0.0_wp
    645    
     628!-- This is necessary because they may not yet have been initialized when they are called from
     629!-- flow_statistics (or - depending on the chosen model run - are never initialized)
     630    sums_divnew_l  = 0.0_wp
     631    sums_divold_l  = 0.0_wp
     632    sums_l_l       = 0.0_wp
     633    sums_wsts_bc_l = 0.0_wp
     634
    646635!
    647636!-- Initialize model variables
    648     IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
     637    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                                &
    649638         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
    650639!
     
    653642          IF ( debug_output )  CALL debug_message( 'initializing with INIFOR', 'start' )
    654643!
    655 !--       Read initial 1D profiles or 3D data from NetCDF file, depending
    656 !--       on the provided level-of-detail.
    657 !--       At the moment, only u, v, w, pt and q are provided. 
     644!--       Read initial 1D profiles or 3D data from NetCDF file, depending on the provided
     645!--       level-of-detail.
     646!--       At the moment, only u, v, w, pt and q are provided.
    658647          CALL netcdf_data_input_init_3d
    659648!
    660 !--       Please note, Inifor provides data from nzb+1 to nzt. 
    661 !--       Bottom and top boundary conditions for Inifor profiles are already
    662 !--       set (just after reading), so that this is not necessary here.
    663 !--       Depending on the provided level-of-detail, initial Inifor data is
    664 !--       either stored on data type (lod=1), or directly on 3D arrays (lod=2).
    665 !--       In order to obtain also initial profiles in case of lod=2 (which
    666 !--       is required for e.g. damping), average over 3D data.
     649!--       Please note, Inifor provides data from nzb+1 to nzt.
     650!--       Bottom and top boundary conditions for Inifor profiles are already set (just after
     651!--       reading), so that this is not necessary here.
     652!--       Depending on the provided level-of-detail, initial Inifor data is either stored on data
     653!--       type (lod=1), or directly on 3D arrays (lod=2).
     654!--       In order to obtain also initial profiles in case of lod=2 (which is required for e.g.
     655!--       damping), average over 3D data.
    667656          IF( init_3d%lod_u == 1 )  THEN
    668657             u_init = init_3d%u_init
    669           ELSEIF( init_3d%lod_u == 2 )  THEN 
    670              ALLOCATE( init_l(nzb:nzt+1) ) 
     658          ELSEIF( init_3d%lod_u == 2 )  THEN
     659             ALLOCATE( init_l(nzb:nzt+1) )
    671660             DO  k = nzb, nzt+1
    672661                init_l(k) = SUM( u(k,nys:nyn,nxl:nxr) )
     
    675664
    676665#if defined( __parallel )
    677              CALL MPI_ALLREDUCE( init_l, u_init, nzt+1-nzb+1,                  &
    678                                  MPI_REAL, MPI_SUM, comm2d, ierr )
     666             CALL MPI_ALLREDUCE( init_l, u_init, nzt+1-nzb+1, MPI_REAL, MPI_SUM, comm2d, ierr )
    679667#else
    680668             u_init = init_l
     
    683671
    684672          ENDIF
    685            
    686           IF( init_3d%lod_v == 1 )  THEN 
     673
     674          IF( init_3d%lod_v == 1 )  THEN
    687675             v_init = init_3d%v_init
    688           ELSEIF( init_3d%lod_v == 2 )  THEN 
    689              ALLOCATE( init_l(nzb:nzt+1) ) 
     676          ELSEIF( init_3d%lod_v == 2 )  THEN
     677             ALLOCATE( init_l(nzb:nzt+1) )
    690678             DO  k = nzb, nzt+1
    691679                init_l(k) = SUM( v(k,nys:nyn,nxl:nxr) )
     
    694682
    695683#if defined( __parallel )
    696              CALL MPI_ALLREDUCE( init_l, v_init, nzt+1-nzb+1,                  &
    697                                  MPI_REAL, MPI_SUM, comm2d, ierr )
     684             CALL MPI_ALLREDUCE( init_l, v_init, nzt+1-nzb+1, MPI_REAL, MPI_SUM, comm2d, ierr )
    698685#else
    699686             v_init = init_l
     
    704691             IF( init_3d%lod_pt == 1 )  THEN
    705692                pt_init = init_3d%pt_init
    706              ELSEIF( init_3d%lod_pt == 2 )  THEN 
    707                 ALLOCATE( init_l(nzb:nzt+1) ) 
     693             ELSEIF( init_3d%lod_pt == 2 )  THEN
     694                ALLOCATE( init_l(nzb:nzt+1) )
    708695                DO  k = nzb, nzt+1
    709696                   init_l(k) = SUM( pt(k,nys:nyn,nxl:nxr) )
     
    712699
    713700#if defined( __parallel )
    714                 CALL MPI_ALLREDUCE( init_l, pt_init, nzt+1-nzb+1,               &
    715                                     MPI_REAL, MPI_SUM, comm2d, ierr )
     701                CALL MPI_ALLREDUCE( init_l, pt_init, nzt+1-nzb+1, MPI_REAL, MPI_SUM, comm2d, ierr )
    716702#else
    717703                pt_init = init_l
     
    725711             IF( init_3d%lod_q == 1 )  THEN
    726712                q_init = init_3d%q_init
    727              ELSEIF( init_3d%lod_q == 2 )  THEN 
    728                 ALLOCATE( init_l(nzb:nzt+1) ) 
     713             ELSEIF( init_3d%lod_q == 2 )  THEN
     714                ALLOCATE( init_l(nzb:nzt+1) )
    729715                DO  k = nzb, nzt+1
    730716                   init_l(k) = SUM( q(k,nys:nyn,nxl:nxr) )
     
    733719
    734720#if defined( __parallel )
    735                 CALL MPI_ALLREDUCE( init_l, q_init, nzt+1-nzb+1,               &
    736                                     MPI_REAL, MPI_SUM, comm2d, ierr )
     721                CALL MPI_ALLREDUCE( init_l, q_init, nzt+1-nzb+1, MPI_REAL, MPI_SUM, comm2d, ierr )
    737722#else
    738723                q_init = init_l
     
    743728
    744729!
    745 !--       Write initial profiles onto 3D arrays.
    746 !--       Work-around, 3D initialization of u,v,w creates artificial
    747 !--       structures wich correlate with the processor grid. The reason
    748 !--       for this is still unknown. To work-around this, 3D initialization
    749 !--       will be effectively reduce to a 1D initialization where no such
    750 !--       artificial structures appear.
     730!--       Write initial profiles onto 3D arrays.
     731!--       Work-around, 3D initialization of u,v,w creates artificial structures which correlate with
     732!--       the processor grid. The reason for this is still unknown. To work-around this, 3D
     733!--       initialization will be effectively reduce to a 1D initialization where no such artificial
     734!--       structures appear.
    751735          DO  i = nxlg, nxrg
    752736             DO  j = nysg, nyng
    753                 IF( init_3d%lod_u == 1  .OR.  init_3d%lod_u == 2 )             &
    754                    u(:,j,i) = u_init(:)
    755                 IF( init_3d%lod_v == 1  .OR.  init_3d%lod_u == 2 )             &
    756                    v(:,j,i) = v_init(:)
    757                 IF( .NOT. neutral  .AND.                                       &
    758                     ( init_3d%lod_pt == 1  .OR.  init_3d%lod_pt == 2 ) )       &
     737                IF( init_3d%lod_u == 1  .OR.  init_3d%lod_u == 2 )  u(:,j,i) = u_init(:)
     738                IF( init_3d%lod_v == 1  .OR.  init_3d%lod_u == 2 )  v(:,j,i) = v_init(:)
     739                IF( .NOT. neutral  .AND.  ( init_3d%lod_pt == 1  .OR.  init_3d%lod_pt == 2 ) )     &
    759740                   pt(:,j,i) = pt_init(:)
    760                 IF( humidity  .AND.                                            &
    761                     ( init_3d%lod_q == 1  .OR.  init_3d%lod_q == 2 ) )         &
     741                IF( humidity  .AND.  ( init_3d%lod_q == 1  .OR.  init_3d%lod_q == 2 ) )            &
    762742                   q(:,j,i) = q_init(:)
    763743             ENDDO
    764744          ENDDO
    765745!
    766 !--       Set geostrophic wind components. 
     746!--       Set geostrophic wind components.
    767747          IF ( init_3d%from_file_ug )  THEN
    768748             ug(:) = init_3d%ug_init(:)
     
    790770
    791771!
    792 !--       Set velocity components at non-atmospheric / oceanic grid points to
    793 !--       zero.
     772!--       Set velocity components at non-atmospheric / oceanic grid points to zero.
    794773          u = MERGE( u, 0.0_wp, BTEST( wall_flags_total_0, 1 ) )
    795774          v = MERGE( v, 0.0_wp, BTEST( wall_flags_total_0, 2 ) )
    796775          w = MERGE( w, 0.0_wp, BTEST( wall_flags_total_0, 3 ) )
    797776!
    798 !--       Initialize surface variables, e.g. friction velocity, momentum
    799 !--       fluxes, etc.
    800           CALL init_surfaces
    801 
    802           IF ( debug_output )  CALL debug_message( 'initializing with INIFOR', 'end' )
     777!--       Initialize surface variables, e.g. friction velocity, momentum fluxes, etc.
     778          CALL  init_surfaces
     779
     780          IF ( debug_output )  CALL  debug_message( 'initializing with INIFOR', 'end' )
    803781!
    804782!--    Initialization via computed 1D-model profiles
    805783       ELSEIF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
    806784
    807           IF ( debug_output )  CALL debug_message( 'initializing with 1D model profiles', 'start' )
     785          IF ( debug_output )  CALL  debug_message( 'initializing with 1D model profiles', 'start' )
    808786!
    809787!--       Use solutions of the 1D model as initial profiles,
     
    833811                   s(:,j,i) = s_init
    834812                ENDDO
    835              ENDDO   
     813             ENDDO
    836814          ENDIF
    837815!
     
    843821!--       Set velocities back to zero
    844822          u = MERGE( u, 0.0_wp, BTEST( wall_flags_total_0, 1 ) )
    845           v = MERGE( v, 0.0_wp, BTEST( wall_flags_total_0, 2 ) )         
    846 !
    847 !--       WARNING: The extra boundary conditions set after running the
    848 !--       -------  1D model impose an error on the divergence one layer
    849 !--                below the topography; need to correct later
    850 !--       ATTENTION: Provisional correction for Piacsek & Williams
    851 !--       ---------  advection scheme: keep u and v zero one layer below
    852 !--                  the topography.
     823          v = MERGE( v, 0.0_wp, BTEST( wall_flags_total_0, 2 ) )
     824!
     825!--       WARNING: The extra boundary conditions set after running the 1D model impose an error on
     826!--       -------- the divergence one layer below the topography; need to correct later
     827!--       ATTENTION: Provisional correction for Piacsek & Williams advection scheme: keep u and v
     828!--       ---------- zero one layer below the topography.
    853829          IF ( ibc_uv_b == 1 )  THEN
    854830!
     
    863839          ENDIF
    864840!
    865 !--       Initialize surface variables, e.g. friction velocity, momentum
    866 !--       fluxes, etc.
     841!--       Initialize surface variables, e.g. friction velocity, momentum fluxes, etc.
    867842          CALL init_surfaces
    868843
    869           IF ( debug_output )  CALL debug_message( 'initializing with 1D model profiles', 'end' )
    870 
    871        ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
    872        THEN
    873 
    874           IF ( debug_output )  CALL debug_message( 'initializing with constant profiles', 'start' )
    875 
    876 !
    877 !--       Use constructed initial profiles (velocity constant with height,
    878 !--       temperature profile with constant gradient)
     844          IF ( debug_output )  CALL  debug_message( 'initializing with 1D model profiles', 'end' )
     845
     846       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )  THEN
     847
     848          IF ( debug_output )  CALL  debug_message( 'initializing with constant profiles', 'start' )
     849
     850!
     851!--       Use constructed initial profiles (velocity constant with height, temperature profile with
     852!--       constant gradient)
    879853          DO  i = nxlg, nxrg
    880854             DO  j = nysg, nyng
     
    889863          v = MERGE( v, 0.0_wp, BTEST( wall_flags_total_0, 2 ) )
    890864!
    891 !--       Set initial horizontal velocities at the lowest computational grid
    892 !--       levels to zero in order to avoid too small time steps caused by the
    893 !--       diffusion limit in the initial phase of a run (at k=1, dz/2 occurs
    894 !--       in the limiting formula!).
    895 !--       Please note, in case land- or urban-surface model is used and a
    896 !--       spinup is applied, masking the lowest computational level is not
    897 !--       possible as MOST as well as energy-balance parametrizations will not
    898 !--       work with zero wind velocity.
    899           IF ( ibc_uv_b /= 1  .AND.  .NOT.  spinup )  THEN
     865!--       Set initial horizontal velocities at the lowest computational grid levels to zero in order
     866!--       to avoid too small time steps caused by the diffusion limit in the initial phase of a run
     867!--       (at k=1, dz/2 occurs in the limiting formula!).
     868!--       Please note, in case land- or urban-surface model is used and a spinup is applied, masking
     869!--       the lowest computational level is not possible as MOST as well as energy-balance
     870!--       parametrizations will not work with zero wind velocity.
     871          IF ( ibc_uv_b /= 1  .AND.  .NOT. spinup )  THEN
    900872             DO  i = nxlg, nxrg
    901873                DO  j = nysg, nyng
    902874                   DO  k = nzb, nzt
    903                       u(k,j,i) = MERGE( u(k,j,i), 0.0_wp,                      &
    904                                      BTEST( wall_flags_total_0(k,j,i), 20 ) )
    905                       v(k,j,i) = MERGE( v(k,j,i), 0.0_wp,                      &
    906                                      BTEST( wall_flags_total_0(k,j,i), 21 ) )
     875                      u(k,j,i) = MERGE( u(k,j,i), 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 20 ) )
     876                      v(k,j,i) = MERGE( v(k,j,i), 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 21 ) )
    907877                   ENDDO
    908878                ENDDO
     
    917887             ENDDO
    918888          ENDIF
    919          
     889
    920890          IF ( passive_scalar )  THEN
    921891             DO  i = nxlg, nxrg
     
    927897
    928898!
    929 !--       Compute initial temperature field and other constants used in case
    930 !--       of a sloping surface
     899!--       Compute initial temperature field and other constants used in case of a sloping surface.
    931900          IF ( sloping_surface )  CALL init_slope
    932901!
    933 !--       Initialize surface variables, e.g. friction velocity, momentum
    934 !--       fluxes, etc.
     902!--       Initialize surface variables, e.g. friction velocity, momentum fluxes, etc.
    935903          CALL init_surfaces
    936          
     904
    937905          IF ( debug_output )  CALL debug_message( 'initializing with constant profiles', 'end' )
    938906
    939        ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 )                  &
    940        THEN
     907       ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 )  THEN
    941908
    942909          IF ( debug_output )  CALL debug_message( 'initializing by user', 'start' )
    943910!
    944 !--       Pre-initialize surface variables, i.e. setting start- and end-indices
    945 !--       at each (j,i)-location. Please note, this does not supersede
    946 !--       user-defined initialization of surface quantities.
     911!--       Pre-initialize surface variables, i.e. setting start- and end-indices at each
     912!--       (j,i)-location. Please note, this does not supersede user-defined initialization of
     913!--       surface quantities.
    947914          CALL init_surfaces
    948915!
     
    954921       ENDIF
    955922
    956        IF ( debug_output )  CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'start' )
     923       IF ( debug_output )  THEN
     924          CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'start' )
     925       ENDIF
    957926
    958927!
    959928!--    Bottom boundary
    960        IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2  )  THEN
     929       IF ( ibc_uv_b == 0  .OR. ibc_uv_b == 2  )  THEN
    961930          u(nzb,:,:) = 0.0_wp
    962931          v(nzb,:,:) = 0.0_wp
     
    975944
    976945!
    977 !--    Store initial profiles for output purposes etc.. Please note, in case of
    978 !--    initialization of u, v, w, pt, and q via output data derived from larger
    979 !--    scale models, data will not be horizontally homogeneous. Actually, a mean
    980 !--    profile should be calculated before.   
     946!--    Store initial profiles for output purposes etc.. Please note, in case of initialization of u,
     947!--    v, w, pt, and q via output data derived from larger scale models, data will not be
     948!--    horizontally homogeneous. Actually, a mean profile should be calculated before.
    981949       hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
    982950       hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
     
    989957       IF ( humidity )  THEN
    990958!
    991 !--       Store initial profile of total water content, virtual potential
    992 !--       temperature
     959!--       Store initial profile of total water content, virtual potential temperature
    993960          hom(:,1,26,:) = SPREAD(   q(:,nys,nxl), 2, statistic_regions+1 )
    994961          hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
    995962!
    996 !--       Store initial profile of mixing ratio and potential
    997 !--       temperature
     963!--       Store initial profile of mixing ratio and potential temperature
    998964          IF ( bulk_cloud_model  .OR.  cloud_droplets ) THEN
    999965             hom(:,1,27,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
     
    1011977!--    Initialize the random number generators (from numerical recipes)
    1012978       CALL random_function_ini
    1013        
     979
    1014980       IF ( random_generator == 'random-parallel' )  THEN
    1015981          CALL init_parallel_random_generator( nx, nys, nyn, nxl, nxr )
    1016982       ENDIF
    1017983!
    1018 !--    Set the reference state to be used in the buoyancy terms (for ocean runs
    1019 !--    the reference state will be set (overwritten) in init_ocean)
     984!--    Set the reference state to be used in the buoyancy terms (for ocean runs the reference state
     985!--    will be set (overwritten) in init_ocean).
    1020986       IF ( use_single_reference_value )  THEN
    1021           IF (  .NOT. humidity )  THEN
     987          IF ( .NOT. humidity )  THEN
    1022988             ref_state(:) = pt_reference
    1023989          ELSE
     
    1025991          ENDIF
    1026992       ELSE
    1027           IF (  .NOT. humidity )  THEN
     993          IF ( .NOT. humidity )  THEN
    1028994             ref_state(:) = pt_init(:)
    1029995          ELSE
     
    10501016
    10511017!
    1052 !--    Impose temperature anomaly (advection test only) or warm air bubble
    1053 !--    close to surface
    1054        IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0  .OR.  &
     1018!--    Impose temperature anomaly (advection test only) or warm air bubble close to surface.
     1019       IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0  .OR.                          &
    10551020            INDEX( initializing_actions, 'initialize_bubble' ) /= 0  )  THEN
    10561021          CALL init_pt_anomaly
    10571022       ENDIF
    1058        
     1023
    10591024!
    10601025!--    If required, change the surface temperature at the start of the 3D run
     
    10661031!--    If required, change the surface humidity/scalar at the start of the 3D
    10671032!--    run
    1068        IF ( humidity  .AND.  q_surface_initial_change /= 0.0_wp )              &
     1033       IF ( humidity  .AND.  q_surface_initial_change /= 0.0_wp )                                  &
    10691034          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
    1070          
    1071        IF ( passive_scalar .AND.  s_surface_initial_change /= 0.0_wp )         &
     1035
     1036       IF ( passive_scalar  .AND.  s_surface_initial_change /= 0.0_wp )                            &
    10721037          s(nzb,:,:) = s(nzb,:,:) + s_surface_initial_change
    1073        
     1038
    10741039
    10751040!
     
    10821047          q_p = q
    10831048       ENDIF
    1084        
     1049
    10851050       IF ( passive_scalar )  THEN
    10861051          ts_m = 0.0_wp
    10871052          s_p  = s
    1088        ENDIF       
    1089 
    1090        IF ( debug_output )  CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'end' )
     1053       ENDIF
     1054
     1055       IF ( debug_output )  THEN
     1056          CALL debug_message( 'initializing statistics, boundary conditions, etc.', 'end' )
     1057       ENDIF
    10911058
    10921059    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.         &
     
    10941061    THEN
    10951062
    1096        IF ( debug_output )  CALL debug_message( 'initializing in case of restart / cyclic_fill', 'start' )
    1097 !
    1098 !--    Initialize surface elements and its attributes, e.g. heat- and
    1099 !--    momentumfluxes, roughness, scaling parameters. As number of surface
    1100 !--    elements might be different between runs, e.g. in case of cyclic fill,
    1101 !--    and not all surface elements are read, surface elements need to be
    1102 !--    initialized before.
    1103 !--    Please note, in case of cyclic fill, surfaces should be initialized
    1104 !--    after restart data is read, else, individual settings of surface
    1105 !--    parameters will be overwritten from data of precursor run, hence,
    1106 !--    init_surfaces is called a second time after reading the restart data.
    1107        CALL init_surfaces                       
    1108 !
    1109 !--    When reading data for cyclic fill of 3D prerun data files, read
    1110 !--    some of the global variables from the restart file which are required
    1111 !--    for initializing the inflow
     1063       IF ( debug_output )  THEN
     1064          CALL debug_message( 'initializing in case of restart / cyclic_fill', 'start' )
     1065       ENDIF
     1066!
     1067!--    Initialize surface elements and its attributes, e.g. heat- and momentumfluxes, roughness,
     1068!--    scaling parameters. As number of surface elements might be different between runs, e.g. in
     1069!--    case of cyclic fill, and not all surface elements are read, surface elements need to be
     1070!--    initialized before.
     1071!--    Please note, in case of cyclic fill, surfaces should be initialized after restart data is
     1072!--    read, else, individual settings of surface parameters will be overwritten from data of
     1073!--    precursor run, hence, init_surfaces is called a second time after reading the restart data.
     1074       CALL init_surfaces
     1075!
     1076!--    When reading data for cyclic fill of 3D prerun data files, read some of the global variables
     1077!--    from the restart file which are required for initializing the inflow
    11121078       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
    11131079
     
    11331099#endif
    11341100       ENDDO
    1135        
    1136        
     1101
     1102
    11371103       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
    11381104
    11391105!
    1140 !--       In case of cyclic fill, call init_surfaces a second time, so that
    1141 !--       surface properties such as heat fluxes are initialized as prescribed.
     1106!--       In case of cyclic fill, call init_surfaces a second time, so that surface properties such
     1107!--       as heat fluxes are initialized as prescribed.
    11421108          CALL init_surfaces
    11431109
    11441110!
    1145 !--       Overwrite u_init, v_init, pt_init, q_init and s_init with the
    1146 !--       horizontally mean (hom) vertical profiles from the end
    1147 !--       of the prerun, because these profiles shall be used as the reference
    1148 !--       state for the rayleigh damping and the pt_damping. This is especially
    1149 !--       important for the use of large_scale_subsidence, because the
    1150 !--       reference temperature in the free atmosphere changes in time.
     1111!--       Overwrite u_init, v_init, pt_init, q_init and s_init with the horizontally mean (hom)
     1112!--       vertical profiles from the end of the prerun, because these profiles shall be used as the
     1113!--       reference state for the rayleigh damping and the pt_damping. This is especially important
     1114!--       for the use of large_scale_subsidence, because the reference temperature in the free
     1115!--       atmosphere changes in time.
    11511116          u_init(:) = hom_sum(:,1,0)
    11521117          v_init(:) = hom_sum(:,2,0)
    11531118          pt_init(:) = hom_sum(:,4,0)
    1154           IF ( humidity )                                                      &
    1155              q_init(:) = hom_sum(:,41,0)
    1156           IF ( passive_scalar )                                                &
    1157              s_init(:) = hom_sum(:,115,0)
    1158        ENDIF
    1159 !
    1160 !--    In case of complex terrain and cyclic fill method as initialization,
    1161 !--    shift initial data in the vertical direction for each point in the
    1162 !--    x-y-plane depending on local surface height
    1163        IF ( complex_terrain  .AND.                                             &
    1164             TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
     1119          IF ( humidity )  q_init(:) = hom_sum(:,41,0)
     1120          IF ( passive_scalar )  s_init(:) = hom_sum(:,115,0)
     1121       ENDIF
     1122!
     1123!--    In case of complex terrain and cyclic fill method as initialization, shift initial data in
     1124!--    the vertical direction for each point in the x-y-plane depending on local surface height.
     1125       IF ( complex_terrain  .AND.  TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
    11651126          DO  i = nxlg, nxrg
    11661127             DO  j = nysg, nyng
     
    11701131                nz_s_shift = topo_top_ind(j,i,0)
    11711132
    1172                 u(nz_u_shift:nzt+1,j,i)  = u(0:nzt+1-nz_u_shift,j,i)               
     1133                u(nz_u_shift:nzt+1,j,i)  = u(0:nzt+1-nz_u_shift,j,i)
    11731134
    11741135                v(nz_v_shift:nzt+1,j,i)  = v(0:nzt+1-nz_v_shift,j,i)
     
    11831144!
    11841145!--    Initialization of the turbulence recycling method
    1185        IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.               &
    1186             turbulent_inflow )  THEN
     1146       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.  turbulent_inflow )  THEN
    11871147!
    11881148!--       First store the profiles to be used at the inflow.
    1189 !--       These profiles are the (temporally) and horizontally averaged vertical
    1190 !--       profiles from the prerun. Alternatively, prescribed profiles
    1191 !--       for u,v-components can be used.
     1149!--       These profiles are the (temporally) and horizontally averaged vertical profiles from the
     1150!--       prerun. Alternatively, prescribed profiles for u,v-components can be used.
    11921151          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
    11931152
     
    12001159          ENDIF
    12011160          mean_inflow_profiles(:,4) = hom_sum(:,4,0)       ! pt
    1202           IF ( humidity )                                                      &
    1203              mean_inflow_profiles(:,6) = hom_sum(:,41,0)   ! q
    1204           IF ( passive_scalar )                                                &
    1205              mean_inflow_profiles(:,7) = hom_sum(:,115,0)   ! s
    1206 
    1207 !
    1208 !--       In case of complex terrain, determine vertical displacement at inflow
    1209 !--       boundary and adjust mean inflow profiles
     1161          IF ( humidity )  mean_inflow_profiles(:,6) = hom_sum(:,41,0)          ! q
     1162          IF ( passive_scalar )  mean_inflow_profiles(:,7) = hom_sum(:,115,0)   ! s
     1163
     1164!
     1165!--       In case of complex terrain, determine vertical displacement at inflow boundary and adjust
     1166!--       mean inflow profiles
    12101167          IF ( complex_terrain )  THEN
    1211              IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 )  THEN
     1168             IF ( nxlg <= 0  .AND.  nxrg >= 0  .AND.  nysg <= 0  .AND. nyng >= 0 )  THEN
    12121169                nz_u_shift_l = topo_top_ind(j,i,1)
    12131170                nz_v_shift_l = topo_top_ind(j,i,2)
     
    12221179
    12231180#if defined( __parallel )
    1224              CALL MPI_ALLREDUCE(nz_u_shift_l, nz_u_shift, 1, MPI_INTEGER,      &
    1225                                 MPI_MAX, comm2d, ierr)
    1226              CALL MPI_ALLREDUCE(nz_v_shift_l, nz_v_shift, 1, MPI_INTEGER,      &
    1227                                 MPI_MAX, comm2d, ierr)
    1228              CALL MPI_ALLREDUCE(nz_w_shift_l, nz_w_shift, 1, MPI_INTEGER,      &
    1229                                 MPI_MAX, comm2d, ierr)
    1230              CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER,      &
    1231                                 MPI_MAX, comm2d, ierr)
     1181             CALL MPI_ALLREDUCE( nz_u_shift_l, nz_u_shift, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
     1182             CALL MPI_ALLREDUCE( nz_v_shift_l, nz_v_shift, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
     1183             CALL MPI_ALLREDUCE( nz_w_shift_l, nz_w_shift, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
     1184             CALL MPI_ALLREDUCE( nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
    12321185#else
    12331186             nz_u_shift = nz_u_shift_l
     
    12481201
    12491202!
    1250 !--       If necessary, adjust the horizontal flow field to the prescribed
    1251 !--       profiles
     1203!--       If necessary, adjust the horizontal flow field to the prescribed profiles
    12521204          IF ( use_prescribed_profile_data )  THEN
    12531205             DO  i = nxlg, nxrg
     
    12621214
    12631215!
    1264 !--       Use these mean profiles at the inflow (provided that Dirichlet
    1265 !--       conditions are used)
     1216!--       Use these mean profiles at the inflow (provided that Dirichlet conditions are used)
    12661217          IF ( bc_dirichlet_l )  THEN
    12671218             DO  j = nysg, nyng
     
    12711222                   w(k,j,nxlg:-1)  = 0.0_wp
    12721223                   pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
    1273                    IF ( humidity )                                             &
    1274                       q(k,j,nxlg:-1)  = mean_inflow_profiles(k,6)
    1275                    IF ( passive_scalar )                                       &
    1276                       s(k,j,nxlg:-1)  = mean_inflow_profiles(k,7)                     
     1224                   IF ( humidity )  q(k,j,nxlg:-1)  = mean_inflow_profiles(k,6)
     1225                   IF ( passive_scalar )  s(k,j,nxlg:-1)  = mean_inflow_profiles(k,7)
    12771226                ENDDO
    12781227             ENDDO
     
    12801229
    12811230!
    1282 !--       Calculate the damping factors to be used at the inflow. For a
    1283 !--       turbulent inflow the turbulent fluctuations have to be limited
    1284 !--       vertically because otherwise the turbulent inflow layer will grow
    1285 !--       in time.
     1231!--       Calculate the damping factors to be used at the inflow. For a turbulent inflow the
     1232!--       turbulent fluctuations have to be limited vertically because otherwise the turbulent
     1233!--       inflow layer will grow in time.
    12861234          IF ( inflow_damping_height == 9999999.9_wp )  THEN
    12871235!
    1288 !--          Default: use the inversion height calculated by the prerun; if
    1289 !--          this is zero, inflow_damping_height must be explicitly
    1290 !--          specified.
     1236!--          Default: use the inversion height calculated by the prerun; if this is zero,
     1237!--          inflow_damping_height must be explicitly specified.
    12911238             IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0_wp )  THEN
    12921239                inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
    12931240             ELSE
    1294                 WRITE( message_string, * ) 'inflow_damping_height must be ',   &
    1295                      'explicitly specified because&the inversion height ',     &
    1296                      'calculated by the prerun is zero.'
     1241                WRITE( message_string, * ) 'inflow_damping_height must be ',                       &
     1242                                           'explicitly specified because&the inversion height ',   &
     1243                                           'calculated by the prerun is zero.'
    12971244                CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
    12981245             ENDIF
     
    13021249          IF ( inflow_damping_width == 9999999.9_wp )  THEN
    13031250!
    1304 !--          Default for the transition range: one tenth of the undamped
    1305 !--          layer
     1251!--          Default for the transition range: one tenth of the undamped layer
    13061252             inflow_damping_width = 0.1_wp * inflow_damping_height
    13071253
     
    13151261                inflow_damping_factor(k) = 1.0_wp
    13161262             ELSEIF ( zu(k) <= ( inflow_damping_height + inflow_damping_width ) )  THEN
    1317                 inflow_damping_factor(k) = 1.0_wp -                            &
    1318                                            ( zu(k) - inflow_damping_height ) / &
    1319                                            inflow_damping_width
     1263                inflow_damping_factor(k) = 1.0_wp -                                                &
     1264                                           ( zu(k) - inflow_damping_height ) / inflow_damping_width
    13201265             ELSE
    13211266                inflow_damping_factor(k) = 0.0_wp
     
    13281273!
    13291274!--    Inside buildings set velocities back to zero
    1330        IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.                &
    1331             topography /= 'flat' )  THEN
     1275       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.  topography /= 'flat' )  THEN
    13321276!
    13331277!--       Inside buildings set velocities back to zero.
     
    13371281             DO  j = nysg, nyng
    13381282                DO  k = nzb, nzt
    1339                    u(k,j,i)     = MERGE( u(k,j,i), 0.0_wp,                     &
    1340                                       BTEST( wall_flags_total_0(k,j,i), 1 ) )
    1341                    v(k,j,i)     = MERGE( v(k,j,i), 0.0_wp,                     &
    1342                                       BTEST( wall_flags_total_0(k,j,i), 2 ) )
    1343                    w(k,j,i)     = MERGE( w(k,j,i), 0.0_wp,                     &
    1344                                       BTEST( wall_flags_total_0(k,j,i), 3 ) )
     1283                   u(k,j,i) = MERGE( u(k,j,i), 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
     1284                   v(k,j,i) = MERGE( v(k,j,i), 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
     1285                   w(k,j,i) = MERGE( w(k,j,i), 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
    13451286                ENDDO
    13461287             ENDDO
     
    13501291
    13511292!
    1352 !--    Calculate initial temperature field and other constants used in case
    1353 !--    of a sloping surface
     1293!--    Calculate initial temperature field and other constants used in case of a sloping surface
    13541294       IF ( sloping_surface )  CALL init_slope
    13551295
    13561296!
    1357 !--    Initialize new time levels (only done in order to set boundary values
    1358 !--    including ghost points)
     1297!--    Initialize new time levels (only done in order to set boundary values including ghost points)
    13591298       pt_p = pt; u_p = u; v_p = v; w_p = w
    13601299       IF ( humidity )  THEN
     
    13631302       IF ( passive_scalar )  s_p  = s
    13641303!
    1365 !--    Allthough tendency arrays are set in prognostic_equations, they have
    1366 !--    have to be predefined here because they are used (but multiplied with 0)
    1367 !--    there before they are set.
     1304!--    Allthough tendency arrays are set in prognostic_equations, they have have to be predefined
     1305!--    here because they are used (but multiplied with 0) there before they are set.
    13681306       tpt_m = 0.0_wp; tu_m = 0.0_wp; tv_m = 0.0_wp; tw_m = 0.0_wp
    13691307       IF ( humidity )  THEN
     
    13721310       IF ( passive_scalar )  ts_m  = 0.0_wp
    13731311
    1374        IF ( debug_output )  CALL debug_message( 'initializing in case of restart / cyclic_fill', 'end' )
     1312       IF ( debug_output )  THEN
     1313          CALL debug_message( 'initializing in case of restart / cyclic_fill', 'end' )
     1314       ENDIF
    13751315
    13761316    ELSE
     
    14051345          w_m_n(:,:,:) = w(:,ny-1:ny,:)
    14061346       ENDIF
    1407        
     1347
    14081348    ENDIF
    14091349
     
    14201360             DO  j = nys, nyn
    14211361                DO  k = nzb+1, nzt
    1422                    volume_flow_initial_l(1) = volume_flow_initial_l(1) +         &
    1423                                               u_init(k) * dzw(k)                 &
    1424                                      * MERGE( 1.0_wp, 0.0_wp,                    &
    1425                                           BTEST( wall_flags_total_0(k,j,nxr), 1 )&
    1426                                             )
    1427 
    1428                    volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)     &
    1429                                      * MERGE( 1.0_wp, 0.0_wp,                    &
    1430                                           BTEST( wall_flags_total_0(k,j,nxr), 1 )&
    1431                                             )
     1362                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +                           &
     1363                                              u_init(k) * dzw(k)                                   &
     1364                                              * MERGE( 1.0_wp, 0.0_wp,                             &
     1365                                                       BTEST( wall_flags_total_0(k,j,nxr), 1 )     &
     1366                                                     )
     1367
     1368                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)                       &
     1369                                              * MERGE( 1.0_wp, 0.0_wp,                             &
     1370                                                       BTEST( wall_flags_total_0(k,j,nxr), 1 )     &
     1371                                                     )
    14321372                ENDDO
    14331373             ENDDO
    14341374          ENDIF
    1435          
     1375
    14361376          IF ( nyn == ny )  THEN
    14371377             DO  i = nxl, nxr
    14381378                DO  k = nzb+1, nzt
    1439                    volume_flow_initial_l(2) = volume_flow_initial_l(2) +         &
    1440                                               v_init(k) * dzw(k)                 &       
    1441                                      * MERGE( 1.0_wp, 0.0_wp,                    &
    1442                                           BTEST( wall_flags_total_0(k,nyn,i), 2 )&
    1443                                             )
    1444                    volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)     &       
    1445                                      * MERGE( 1.0_wp, 0.0_wp,                    &
    1446                                           BTEST( wall_flags_total_0(k,nyn,i), 2 )&
    1447                                             )
     1379                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +                           &
     1380                                              v_init(k) * dzw(k)                                   &
     1381                                              * MERGE( 1.0_wp, 0.0_wp,                             &
     1382                                                       BTEST( wall_flags_total_0(k,nyn,i), 2 )     &
     1383                                                     )
     1384                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)                       &
     1385                                              * MERGE( 1.0_wp, 0.0_wp,                             &
     1386                                                       BTEST( wall_flags_total_0(k,nyn,i), 2 )     &
     1387                                                     )
    14481388                ENDDO
    14491389             ENDDO
     
    14511391
    14521392#if defined( __parallel )
    1453           CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
    1454                               2, MPI_REAL, MPI_SUM, comm2d, ierr )
    1455           CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
    1456                               2, MPI_REAL, MPI_SUM, comm2d, ierr )
     1393          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1), 2, MPI_REAL,       &
     1394                              MPI_SUM, comm2d, ierr )
     1395          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), 2, MPI_REAL, MPI_SUM,    &
     1396                              comm2d, ierr )
    14571397
    14581398#else
    14591399          volume_flow_initial = volume_flow_initial_l
    14601400          volume_flow_area    = volume_flow_area_l
    1461 #endif 
     1401#endif
    14621402
    14631403       ELSEIF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
     
    14691409             DO  j = nys, nyn
    14701410                DO  k = nzb+1, nzt
    1471                    volume_flow_initial_l(1) = volume_flow_initial_l(1) +         &
    1472                                               hom_sum(k,1,0) * dzw(k)            &
    1473                                      * MERGE( 1.0_wp, 0.0_wp,                    &
    1474                                           BTEST( wall_flags_total_0(k,j,nx), 1 ) &
    1475                                             )
    1476                    volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)     &
    1477                                      * MERGE( 1.0_wp, 0.0_wp,                    &
    1478                                           BTEST( wall_flags_total_0(k,j,nx), 1 ) &
    1479                                             )
     1411                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +                           &
     1412                                              hom_sum(k,1,0) * dzw(k)                              &
     1413                                              * MERGE( 1.0_wp, 0.0_wp,                             &
     1414                                                       BTEST( wall_flags_total_0(k,j,nx), 1 )      &
     1415                                                     )
     1416                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)                       &
     1417                                              * MERGE( 1.0_wp, 0.0_wp,                             &
     1418                                                       BTEST( wall_flags_total_0(k,j,nx), 1 )      &
     1419                                                     )
    14801420                ENDDO
    14811421             ENDDO
    14821422          ENDIF
    1483          
     1423
    14841424          IF ( nyn == ny )  THEN
    14851425             DO  i = nxl, nxr
    14861426                DO  k = nzb+1, nzt
    1487                    volume_flow_initial_l(2) = volume_flow_initial_l(2) +         &
    1488                                               hom_sum(k,2,0) * dzw(k)            &       
    1489                                      * MERGE( 1.0_wp, 0.0_wp,                    &
    1490                                           BTEST( wall_flags_total_0(k,ny,i), 2 ) &
    1491                                             )
    1492                    volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)     &       
    1493                                      * MERGE( 1.0_wp, 0.0_wp,                    &
    1494                                           BTEST( wall_flags_total_0(k,ny,i), 2 ) &
    1495                                             )
     1427                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +                           &
     1428                                              hom_sum(k,2,0) * dzw(k)                              &
     1429                                              * MERGE( 1.0_wp, 0.0_wp,                             &
     1430                                                       BTEST( wall_flags_total_0(k,ny,i), 2 )      &
     1431                                                     )
     1432                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)                       &
     1433                                              * MERGE( 1.0_wp, 0.0_wp,                             &
     1434                                                       BTEST( wall_flags_total_0(k,ny,i), 2 )      &
     1435                                                     )
    14961436                ENDDO
    14971437             ENDDO
     
    14991439
    15001440#if defined( __parallel )
    1501           CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
    1502                               2, MPI_REAL, MPI_SUM, comm2d, ierr )
    1503           CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
    1504                               2, MPI_REAL, MPI_SUM, comm2d, ierr )
     1441          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1), 2, MPI_REAL,       &
     1442                              MPI_SUM, comm2d, ierr )
     1443          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), 2, MPI_REAL, MPI_SUM,    &
     1444                              comm2d, ierr )
    15051445
    15061446#else
    15071447          volume_flow_initial = volume_flow_initial_l
    15081448          volume_flow_area    = volume_flow_area_l
    1509 #endif 
     1449#endif
    15101450
    15111451       ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
     
    15171457             DO  j = nys, nyn
    15181458                DO  k = nzb+1, nzt
    1519                    volume_flow_initial_l(1) = volume_flow_initial_l(1) +         &
    1520                                               u(k,j,nx) * dzw(k)                 &
    1521                                      * MERGE( 1.0_wp, 0.0_wp,                    &
    1522                                           BTEST( wall_flags_total_0(k,j,nx), 1 ) &
    1523                                             )
    1524                    volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)     &
    1525                                      * MERGE( 1.0_wp, 0.0_wp,                    &
    1526                                           BTEST( wall_flags_total_0(k,j,nx), 1 ) &
    1527                                             )
     1459                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +                           &
     1460                                              u(k,j,nx) * dzw(k)                                   &
     1461                                              * MERGE( 1.0_wp, 0.0_wp,                             &
     1462                                                       BTEST( wall_flags_total_0(k,j,nx), 1 )      &
     1463                                                     )
     1464                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)                       &
     1465                                              * MERGE( 1.0_wp, 0.0_wp,                             &
     1466                                                   BTEST( wall_flags_total_0(k,j,nx), 1 )          &
     1467                                                     )
    15281468                ENDDO
    15291469             ENDDO
    15301470          ENDIF
    1531          
     1471
    15321472          IF ( nyn == ny )  THEN
    15331473             DO  i = nxl, nxr
    15341474                DO  k = nzb+1, nzt
    1535                    volume_flow_initial_l(2) = volume_flow_initial_l(2) +         &
    1536                                               v(k,ny,i) * dzw(k)                 &       
    1537                                      * MERGE( 1.0_wp, 0.0_wp,                    &
    1538                                           BTEST( wall_flags_total_0(k,ny,i), 2 ) &
    1539                                             )
    1540                    volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)     &       
    1541                                      * MERGE( 1.0_wp, 0.0_wp,                    &
    1542                                           BTEST( wall_flags_total_0(k,ny,i), 2 ) &
    1543                                             )
     1475                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +                           &
     1476                                              v(k,ny,i) * dzw(k)                                   &
     1477                                              * MERGE( 1.0_wp, 0.0_wp,                             &
     1478                                                       BTEST( wall_flags_total_0(k,ny,i), 2 )      &
     1479                                                     )
     1480                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)                       &
     1481                                              * MERGE( 1.0_wp, 0.0_wp,                             &
     1482                                                       BTEST( wall_flags_total_0(k,ny,i), 2 )      &
     1483                                                     )
    15441484                ENDDO
    15451485             ENDDO
     
    15471487
    15481488#if defined( __parallel )
    1549           CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
    1550                               2, MPI_REAL, MPI_SUM, comm2d, ierr )
    1551           CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
    1552                               2, MPI_REAL, MPI_SUM, comm2d, ierr )
     1489          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1), 2, MPI_REAL,       &
     1490                              MPI_SUM, comm2d, ierr )
     1491          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), 2, MPI_REAL, MPI_SUM,    &
     1492                              comm2d, ierr )
    15531493
    15541494#else
    15551495          volume_flow_initial = volume_flow_initial_l
    15561496          volume_flow_area    = volume_flow_area_l
    1557 #endif 
    1558 
    1559        ENDIF
    1560 
    1561 !
    1562 !--    In case of 'bulk_velocity' mode, volume_flow_initial is calculated
    1563 !--    from u|v_bulk instead
     1497#endif
     1498
     1499       ENDIF
     1500
     1501!
     1502!--    In case of 'bulk_velocity' mode, volume_flow_initial is calculated from u|v_bulk instead
    15641503       IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
    15651504          volume_flow_initial(1) = u_bulk * volume_flow_area(1)
     
    15691508    ENDIF
    15701509!
    1571 !-- In the following, surface properties can be further initialized with
    1572 !-- input from static driver file.
    1573 !-- At the moment this affects only default surfaces. For example,
    1574 !-- roughness length or sensible / latent heat fluxes can be initialized
    1575 !-- heterogeneously for default surfaces. Therefore, a generic routine
    1576 !-- from netcdf_data_input_mod is called to read a 2D array.
     1510!-- In the following, surface properties can be further initialized with input from static driver
     1511!-- file.
     1512!-- At the moment this affects only default surfaces. For example, roughness length or sensible /
     1513!-- latent heat fluxes can be initialized heterogeneously for default surfaces. Therefore, a generic
     1514!-- routine from netcdf_data_input_mod is called to read a 2D array.
    15771515    IF ( input_pids_static )  THEN
    15781516!
     
    15831521!--    Open the static input file
    15841522#if defined( __netcdf )
    1585        CALL open_read_file( TRIM( input_file_static ) //                    &
    1586                             TRIM( coupling_char ),                          &
    1587                             pids_id )
    1588                            
     1523       CALL open_read_file( TRIM( input_file_static ) //                                           &
     1524                            TRIM( coupling_char ), pids_id )
     1525
    15891526       CALL inquire_num_variables( pids_id, num_var_pids )
    15901527!
     
    15931530       CALL inquire_variable_names( pids_id, vars_pids )
    15941531!
    1595 !--    Input roughness length. 
     1532!--    Input roughness length.
    15961533       IF ( check_existence( vars_pids, 'z0' ) )  THEN
    15971534!
    15981535!--       Read _FillValue attribute
    1599           CALL get_attribute( pids_id, char_fill, tmp_2d%fill,          &
    1600                               .FALSE., 'z0' )                                 
    1601 !                                                                             
    1602 !--       Read variable                                                       
    1603           CALL get_variable( pids_id, 'z0', tmp_2d%var,                 &
    1604                              nxl, nxr, nys, nyn )                             
    1605 !                                                                             
    1606 !--       Initialize roughness length. Note, z0 will be only initialized at   
    1607 !--       default-type surfaces. At natural or urban z0 is implicitly         
    1608 !--       initialized by the respective parameter lists.                       
    1609 !--       Initialize horizontal surface elements.                             
    1610           CALL init_single_surface_properties( surf_def_h(0)%z0,               &
    1611                                                tmp_2d%var,                     &
    1612                                                surf_def_h(0)%ns,               &
    1613                                                tmp_2d%fill,                    &
    1614                                                surf_def_h(0)%i,                &
    1615                                                surf_def_h(0)%j )               
    1616 !                                                                             
    1617 !--       Initialize roughness also at vertical surface elements.             
    1618 !--       Note, the actual 2D input arrays are only defined on the             
    1619 !--       subdomain. Therefore, pass the index arrays with their respective   
    1620 !--       offset values.                                                       
    1621           DO  l = 0, 3                                                         
    1622              CALL init_single_surface_properties(                              &
    1623                                          surf_def_v(l)%z0,                     &
    1624                                          tmp_2d%var,                           &
    1625                                          surf_def_v(l)%ns,                     &
    1626                                          tmp_2d%fill,                          &
    1627                                          surf_def_v(l)%i + surf_def_v(l)%ioff, &
    1628                                          surf_def_v(l)%j + surf_def_v(l)%joff )
     1536          CALL get_attribute( pids_id, char_fill, tmp_2d%fill, .FALSE., 'z0' )
     1537!
     1538!--       Read variable
     1539          CALL get_variable( pids_id, 'z0', tmp_2d%var, nxl, nxr, nys, nyn )
     1540!
     1541!--       Initialize roughness length. Note, z0 will be only initialized at default-type surfaces.
     1542!--       At natural or urban z0 is implicitly initialized by the respective parameter lists.
     1543!--       Initialize horizontal surface elements.
     1544          CALL init_single_surface_properties( surf_def_h(0)%z0, tmp_2d%var, surf_def_h(0)%ns,     &
     1545                                               tmp_2d%fill, surf_def_h(0)%i, surf_def_h(0)%j )
     1546!
     1547!--       Initialize roughness also at vertical surface elements.
     1548!--       Note, the actual 2D input arrays are only defined on the subdomain. Therefore, pass the
     1549!--       index arrays with their respective offset values.
     1550          DO  l = 0, 3
     1551             CALL init_single_surface_properties( surf_def_v(l)%z0, tmp_2d%var, surf_def_v(l)%ns,  &
     1552                                                  tmp_2d%fill, surf_def_v(l)%i+surf_def_v(l)%ioff, &
     1553                                                  surf_def_v(l)%j+surf_def_v(l)%joff )
    16291554          ENDDO
    1630          
    1631        ENDIF
    1632 !
    1633 !--    Input surface sensible heat flux. 
     1555
     1556       ENDIF
     1557!
     1558!--    Input surface sensible heat flux.
    16341559       IF ( check_existence( vars_pids, 'shf' ) )  THEN
    16351560!
    16361561!--       Read _FillValue attribute
    1637           CALL get_attribute( pids_id, char_fill, tmp_2d%fill,                 &
    1638                               .FALSE., 'shf' )
     1562          CALL get_attribute( pids_id, char_fill, tmp_2d%fill, .FALSE., 'shf' )
    16391563!
    16401564!--       Read variable
    1641           CALL get_variable( pids_id, 'shf', tmp_2d%var,                       &
    1642                              nxl, nxr, nys, nyn )
    1643 !
    1644 !--       Initialize heat flux. Note, shf will be only initialized at
    1645 !--       default-type surfaces. At natural or urban shf is implicitly
    1646 !--       initialized by the respective parameter lists.
     1565          CALL get_variable( pids_id, 'shf', tmp_2d%var, nxl, nxr, nys, nyn )
     1566!
     1567!--       Initialize heat flux. Note, shf will be only initialized at default-type surfaces. At
     1568!--       natural or urban shf is implicitly initialized by the respective parameter lists.
    16471569!--       Initialize horizontal surface elements.
    1648           CALL init_single_surface_properties( surf_def_h(0)%shf,              &
    1649                                                tmp_2d%var,                     &
    1650                                                surf_def_h(0)%ns,               &
    1651                                                tmp_2d%fill,                    &
    1652                                                surf_def_h(0)%i,                &
    1653                                                surf_def_h(0)%j )
     1570          CALL init_single_surface_properties( surf_def_h(0)%shf, tmp_2d%var, surf_def_h(0)%ns,    &
     1571                                               tmp_2d%fill, surf_def_h(0)%i, surf_def_h(0)%j )
    16541572!
    16551573!--       Initialize heat flux also at vertical surface elements.
    1656 !--       Note, the actual 2D input arrays are only defined on the
    1657 !--       subdomain. Therefore, pass the index arrays with their respective
    1658 !--       offset values.
     1574!--       Note, the actual 2D input arrays are only defined on the subdomain. Therefore, pass the
     1575!--       index arrays with their respective offset values.
    16591576          DO  l = 0, 3
    1660              CALL init_single_surface_properties(                              &
    1661                                          surf_def_v(l)%shf,                    &
    1662                                          tmp_2d%var,                           &
    1663                                          surf_def_v(l)%ns,                     &
    1664                                          tmp_2d%fill,                          &
    1665                                          surf_def_v(l)%i + surf_def_v(l)%ioff, &
    1666                                          surf_def_v(l)%j + surf_def_v(l)%joff )
     1577             CALL init_single_surface_properties( surf_def_v(l)%shf, tmp_2d%var, surf_def_v(l)%ns, &
     1578                                                  tmp_2d%fill, surf_def_v(l)%i+surf_def_v(l)%ioff, &
     1579                                                  surf_def_v(l)%j+surf_def_v(l)%joff )
    16671580          ENDDO
    16681581
    16691582       ENDIF
    16701583!
    1671 !--    Input surface sensible heat flux. 
     1584!--    Input surface sensible heat flux.
    16721585       IF ( check_existence( vars_pids, 'qsws' ) )  THEN
    16731586!
     
    16801593                             nxl, nxr, nys, nyn )
    16811594!
    1682 !--       Initialize latent heat flux. Note, qsws will be only initialized at
    1683 !--       default-type surfaces. At natural or urban qsws is implicitly
    1684 !--       initialized by the respective parameter lists.
     1595!--       Initialize latent heat flux. Note, qsws will be only initialized at default-type surfaces.
     1596!--       At natural or urban qsws is implicitly initialized by the respective parameter lists.
    16851597!--       Initialize horizontal surface elements.
    1686           CALL init_single_surface_properties( surf_def_h(0)%qsws,             &
    1687                                                tmp_2d%var,                     &
    1688                                                surf_def_h(0)%ns,               &
    1689                                                tmp_2d%fill,                    &
    1690                                                surf_def_h(0)%i,                &
    1691                                                surf_def_h(0)%j )
     1598          CALL init_single_surface_properties( surf_def_h(0)%qsws, tmp_2d%var, surf_def_h(0)%ns,   &
     1599                                               tmp_2d%fill, surf_def_h(0)%i, surf_def_h(0)%j )
    16921600!
    16931601!--       Initialize latent heat flux also at vertical surface elements.
    1694 !--       Note, the actual 2D input arrays are only defined on the
    1695 !--       subdomain. Therefore, pass the index arrays with their respective
    1696 !--       offset values.
     1602!--       Note, the actual 2D input arrays are only defined on the subdomain. Therefore, pass the
     1603!--       index arrays with their respective offset values.
    16971604          DO  l = 0, 3
    1698              CALL init_single_surface_properties(                              &
    1699                                          surf_def_v(l)%qsws,                   &
    1700                                          tmp_2d%var,                           &
    1701                                          surf_def_v(l)%ns,                     &
    1702                                          tmp_2d%fill,                          &
    1703                                          surf_def_v(l)%i + surf_def_v(l)%ioff, &
    1704                                          surf_def_v(l)%j + surf_def_v(l)%joff )
     1605             CALL init_single_surface_properties( surf_def_v(l)%qsws, tmp_2d%var, surf_def_v(l)%ns,&
     1606                                                  tmp_2d%fill, surf_def_v(l)%i+surf_def_v(l)%ioff, &
     1607                                                  surf_def_v(l)%j+surf_def_v(l)%joff )
    17051608          ENDDO
    17061609
    17071610       ENDIF
    17081611!
    1709 !--    Additional variables, can be initialized the 
     1612!--    Additional variables, can be initialized the
    17101613!--    same way.
    17111614
     
    17131616!--    Finally, close the input file and deallocate temporary arrays
    17141617       DEALLOCATE( vars_pids )
    1715        
     1618
    17161619       CALL close_input_file( pids_id )
    17171620#endif
     
    17191622    ENDIF
    17201623!
    1721 !-- Finally, if random_heatflux is set, disturb shf at horizontal
    1722 !-- surfaces. Actually, this should be done in surface_mod, where all other
    1723 !-- initializations of surface quantities are done. However, this
    1724 !-- would create a ring dependency, hence, it is done here. Maybe delete
    1725 !-- disturb_heatflux and tranfer the respective code directly into the
    1726 !-- initialization in surface_mod.         
    1727     IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
     1624!-- Finally, if random_heatflux is set, disturb shf at horizontal surfaces. Actually, this should be
     1625!-- done in surface_mod, where all other initializations of surface quantities are done. However,
     1626!-- this would create a ring dependency, hence, it is done here. Maybe delete disturb_heatflux and
     1627!-- tranfer the respective code directly into the initialization in surface_mod.
     1628    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                                &
    17281629         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
    1729  
    1730        IF ( use_surface_fluxes  .AND.  constant_heatflux  .AND.                &
    1731             random_heatflux )  THEN
     1630
     1631       IF ( use_surface_fluxes  .AND.  constant_heatflux  .AND.  random_heatflux )  THEN
    17321632          IF ( surf_def_h(0)%ns >= 1 )  CALL disturb_heatflux( surf_def_h(0) )
    17331633          IF ( surf_lsm_h%ns    >= 1 )  CALL disturb_heatflux( surf_lsm_h    )
     
    17371637
    17381638!
    1739 !-- Compute total sum of grid points and the mean surface level height for each
    1740 !-- statistic region. These are mainly used for horizontal averaging of
    1741 !-- turbulence statistics.
    1742 !-- ngp_2dh: number of grid points of a horizontal cross section through the
    1743 !--          respective statistic region
     1639!-- Compute total sum of grid points and the mean surface level height for each statistic region.
     1640!-- These are mainly used for horizontal averaging of turbulence statistics.
     1641!-- ngp_2dh: number of grid points of a horizontal cross section through the respective statistic
     1642!--          region
    17441643!-- ngp_3d:  number of grid points of the respective statistic region
    17451644    ngp_2dh_outer_l   = 0
     
    17661665                ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
    17671666!
    1768 !--             Determine mean surface-level height. In case of downward-
    1769 !--             facing walls are present, more than one surface level exist.
    1770 !--             In this case, use the lowest surface-level height.
    1771                 IF ( surf_def_h(0)%start_index(j,i) <=                         &
    1772                      surf_def_h(0)%end_index(j,i) )  THEN
     1667!--             Determine mean surface-level height. In case of downward-facing walls are present,
     1668!--             more than one surface level exist.
     1669!--             In this case, use the lowest surface-level height.
     1670                IF ( surf_def_h(0)%start_index(j,i) <= surf_def_h(0)%end_index(j,i) )  THEN
    17731671                   m = surf_def_h(0)%start_index(j,i)
    17741672                   k = surf_def_h(0)%k(m)
    1775                    mean_surface_level_height_l(sr) =                           &
    1776                                        mean_surface_level_height_l(sr) + zw(k-1)
     1673                   mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) + zw(k-1)
    17771674                ENDIF
    1778                 IF ( surf_lsm_h%start_index(j,i) <=                            &
    1779                      surf_lsm_h%end_index(j,i) )  THEN
     1675                IF ( surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i) )  THEN
    17801676                   m = surf_lsm_h%start_index(j,i)
    17811677                   k = surf_lsm_h%k(m)
    1782                    mean_surface_level_height_l(sr) =                           &
    1783                                        mean_surface_level_height_l(sr) + zw(k-1)
     1678                   mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) + zw(k-1)
    17841679                ENDIF
    1785                 IF ( surf_usm_h%start_index(j,i) <=                            &
    1786                      surf_usm_h%end_index(j,i) )  THEN
     1680                IF ( surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i) )  THEN
    17871681                   m = surf_usm_h%start_index(j,i)
    17881682                   k = surf_usm_h%k(m)
    1789                    mean_surface_level_height_l(sr) =                           &
    1790                                        mean_surface_level_height_l(sr) + zw(k-1)
     1683                   mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) + zw(k-1)
    17911684                ENDIF
    17921685
     
    17961689!
    17971690!--                xy-grid points above topography
    1798                    ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr)     +         &
    1799                              MERGE( 1, 0, BTEST( wall_flags_total_0(k,j,i), 24 ) )
    1800 
    1801                    ngp_2dh_s_inner_l(k,sr) = ngp_2dh_s_inner_l(k,sr) +         &
    1802                              MERGE( 1, 0, BTEST( wall_flags_total_0(k,j,i), 22 ) )
     1691                   ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr)     +                             &
     1692                                           MERGE( 1, 0, BTEST( wall_flags_total_0(k,j,i), 24 ) )
     1693
     1694                   ngp_2dh_s_inner_l(k,sr) = ngp_2dh_s_inner_l(k,sr) +                             &
     1695                                             MERGE( 1, 0, BTEST( wall_flags_total_0(k,j,i), 22 ) )
    18031696
    18041697                ENDDO
     
    18171710#if defined( __parallel )
    18181711    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1819     CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM,    &
     1712    CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
     1713    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     1714    CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr, MPI_INTEGER, MPI_SUM, &
    18201715                        comm2d, ierr )
    18211716    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1822     CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr,   &
    1823                         MPI_INTEGER, MPI_SUM, comm2d, ierr )
     1717    CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0), (nz+2)*sr, MPI_INTEGER,      &
     1718                        MPI_SUM, comm2d, ierr )
    18241719    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1825     CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0),          &
    1826                         (nz+2)*sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
    1827     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1828     CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL,  &
    1829                         MPI_SUM, comm2d, ierr )
     1720    CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL, MPI_SUM, comm2d,     &
     1721                        ierr )
    18301722    ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) )
    18311723    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1832     CALL MPI_ALLREDUCE( mean_surface_level_height_l(0),                        &
    1833                         mean_surface_level_height(0), sr, MPI_REAL,            &
     1724    CALL MPI_ALLREDUCE( mean_surface_level_height_l(0), mean_surface_level_height(0), sr, MPI_REAL,&
    18341725                        MPI_SUM, comm2d, ierr )
    18351726    mean_surface_level_height = mean_surface_level_height / REAL( ngp_2dh )
     
    18421733#endif
    18431734
    1844     ngp_3d = INT ( ngp_2dh, KIND = SELECTED_INT_KIND( 18 ) ) * &
     1735    ngp_3d = INT ( ngp_2dh, KIND = SELECTED_INT_KIND( 18 ) ) *                                     &
    18451736             INT ( (nz + 2 ), KIND = SELECTED_INT_KIND( 18 ) )
    18461737
    18471738!
    1848 !-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics,
    1849 !-- buoyancy, etc. A zero value will occur for cases where all grid points of
    1850 !-- the respective subdomain lie below the surface topography
    1851     ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   )
    1852     ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )),             &
    1853                            ngp_3d_inner(:) )
    1854     ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) )
    1855 
    1856     DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l,       &
    1857                 ngp_3d_inner_l, ngp_3d_inner_tmp )
    1858 
    1859 !
    1860 !-- Initializae 3D offline nesting in COSMO model and read data from
     1739!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics, buoyancy, etc. A
     1740!-- zero value will occur for cases where all grid points of the respective subdomain lie below the
     1741!-- surface topography
     1742    ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   )
     1743    ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )), ngp_3d_inner(:) )
     1744    ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) )
     1745
     1746    DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l,           &
     1747                ngp_3d_inner_tmp )
     1748
     1749!
     1750!-- Initializae 3D offline nesting in COSMO model and read data from
    18611751!-- external NetCDF file.
    18621752    IF ( nesting_offline )  CALL nesting_offl_init
     
    18681758!-- Impose random perturbation on the horizontal velocity field and then
    18691759!-- remove the divergences from the velocity field at the initial stage
    1870     IF ( create_disturbances  .AND.  disturbance_energy_limit /= 0.0_wp  .AND. &
    1871          TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
     1760    IF ( create_disturbances  .AND.  disturbance_energy_limit /= 0.0_wp  .AND.                     &
     1761         TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                                &
    18721762         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
    18731763
    1874        IF ( debug_output )  CALL debug_message( 'creating disturbances + applying pressure solver', 'start' )
     1764       IF ( debug_output )  THEN
     1765          CALL debug_message( 'creating disturbances + applying pressure solver', 'start' )
     1766       ENDIF
    18751767!
    18761768!--    Needed for both disturb_field and pres
     
    19051797!$ACC END DATA
    19061798
    1907        IF ( debug_output )  CALL debug_message( 'creating disturbances + applying pressure solver', 'end' )
     1799       IF ( debug_output )  THEN
     1800          CALL debug_message( 'creating disturbances + applying pressure solver', 'end' )
     1801       ENDIF
    19081802
    19091803    ENDIF
     
    19181812!--    Check temperature in case of too large domain height
    19191813       DO  k = nzb, nzt+1
    1920           IF ( ( pt_surface * exner_function(surface_pressure * 100.0_wp) - g/c_p * zu(k) ) < 0.0_wp )  THEN
    1921              WRITE( message_string, * )  'absolute temperature < 0.0 at zu(', k, &
    1922                                          ') = ', zu(k)
     1814          IF ( ( pt_surface * exner_function( surface_pressure * 100.0_wp ) - g/c_p * zu(k) )      &
     1815                 < 0.0_wp )  THEN
     1816             WRITE( message_string, * )  'absolute temperature < 0.0 at zu(', k, ') = ', zu(k)
    19231817             CALL message( 'init_3d_model', 'PA0142', 1, 2, 0, 6, 0 )
    19241818          ENDIF
     
    19271821!
    19281822!--    Calculate vertical profile of the hydrostatic pressure (hyp)
    1929        hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
    1930        d_exner = exner_function_invers(hyp)
    1931        exner = 1.0_wp / exner_function_invers(hyp)
    1932        hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
     1823       hyp    = barometric_formula( zu, pt_surface * exner_function( surface_pressure * 100.0_wp ),&
     1824                                    surface_pressure * 100.0_wp )
     1825       d_exner = exner_function_invers( hyp )
     1826       exner = 1.0_wp / exner_function_invers( hyp )
     1827       hyrho  = ideal_gas_law_rho_pt( hyp, pt_init )
    19331828!
    19341829!--    Compute reference density
    1935        rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
     1830       rho_surface = ideal_gas_law_rho( surface_pressure * 100.0_wp,                               &
     1831                                        pt_surface * exner_function( surface_pressure * 100.0_wp ) )
    19361832
    19371833    ENDIF
     
    19471843    CALL module_interface_init
    19481844!
    1949 !-- Initialize surface layer (done after LSM as roughness length are required
    1950 !-- for initialization
     1845!-- Initialize surface layer (done after LSM as roughness length are required for initialization
    19511846    IF ( constant_flux_layer )  CALL init_surface_layer_fluxes
    19521847!
     
    19541849    IF ( surface_output )  CALL surface_data_output_init
    19551850!
    1956 !-- Initialize the ws-scheme.   
    1957     IF ( ws_scheme_sca .OR. ws_scheme_mom )  CALL ws_init
     1851!-- Initialize the ws-scheme.
     1852    IF ( ws_scheme_sca  .OR. ws_scheme_mom )  CALL ws_init
    19581853!
    19591854!-- Perform post-initializing checks for all other modules
     
    19611856
    19621857!
    1963 !-- Initialize surface forcing corresponding to large-scale forcing. Therein, 
     1858!-- Initialize surface forcing corresponding to large-scale forcing. Therein,
    19641859!-- initialize heat-fluxes, etc. via datatype. Revise it later!
    1965     IF ( large_scale_forcing .AND. lsf_surf )  THEN
     1860    IF ( large_scale_forcing  .AND. lsf_surf )  THEN
    19661861       IF ( use_surface_fluxes  .AND.  constant_heatflux )  THEN
    1967           CALL ls_forcing_surf ( simulated_time )
    1968        ENDIF
    1969     ENDIF
    1970 !
    1971 !-- Setting weighting factors for calculation of perturbation pressure
    1972 !-- and turbulent quantities from the RK substeps
    1973     IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN      ! for RK3-method
     1862          CALL ls_forcing_surf( simulated_time )
     1863       ENDIF
     1864    ENDIF
     1865!
     1866!-- Setting weighting factors for calculation of perturbation pressure and turbulent quantities from
     1867!-- the RK substeps.
     1868    IF ( TRIM( timestep_scheme ) == 'runge-kutta-3' )  THEN      ! for RK3-method
    19741869
    19751870       weight_substep(1) = 1._wp/6._wp
     
    19811876       weight_pres(3)    = 1._wp/4._wp
    19821877
    1983     ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! for RK2-method
     1878    ELSEIF ( TRIM( timestep_scheme ) == 'runge-kutta-2' )  THEN  ! for RK2-method
    19841879
    19851880       weight_substep(1) = 1._wp/2._wp
    19861881       weight_substep(2) = 1._wp/2._wp
    1987          
     1882
    19881883       weight_pres(1)    = 1._wp/2._wp
    1989        weight_pres(2)    = 1._wp/2._wp       
     1884       weight_pres(2)    = 1._wp/2._wp
    19901885
    19911886    ELSE                                     ! for Euler-method
    19921887
    1993        weight_substep(1) = 1.0_wp     
    1994        weight_pres(1)    = 1.0_wp                   
     1888       weight_substep(1) = 1.0_wp
     1889       weight_pres(1)    = 1.0_wp
    19951890
    19961891    ENDIF
     
    20051900          DO  k = nzb+1, nzt
    20061901             IF ( zu(k) >= rayleigh_damping_height )  THEN
    2007                 rdf(k) = rayleigh_damping_factor *                             &
    2008                       ( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height ) &
    2009                              / ( zu(nzt) - rayleigh_damping_height ) )         &
    2010                       )**2
     1902                rdf(k) = rayleigh_damping_factor *                                                 &
     1903                         ( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height )                  &
     1904                                / ( zu(nzt) - rayleigh_damping_height ) )                          &
     1905                         )**2
    20111906             ENDIF
    20121907          ENDDO
    20131908       ELSE
    20141909!
    2015 !--       In ocean mode, rayleigh damping is applied in the lower part of the
    2016 !--       model domain
     1910!--       In ocean mode, rayleigh damping is applied in the lower part of the model domain
    20171911          DO  k = nzt, nzb+1, -1
    20181912             IF ( zu(k) <= rayleigh_damping_height )  THEN
    2019                 rdf(k) = rayleigh_damping_factor *                             &
    2020                       ( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) ) &
    2021                              / ( rayleigh_damping_height - zu(nzb+1) ) )       &
    2022                       )**2
     1913                rdf(k) = rayleigh_damping_factor *                                                 &
     1914                         ( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) )                  &
     1915                                / ( rayleigh_damping_height - zu(nzb+1) ) )                        &
     1916                         )**2
    20231917             ENDIF
    20241918          ENDDO
     
    20291923
    20301924!
    2031 !-- Initialize the starting level and the vertical smoothing factor used for
    2032 !-- the external pressure gradient
     1925!-- Initialize the starting level and the vertical smoothing factor used for the external pressure
     1926!-- gradient
    20331927    dp_smooth_factor = 1.0_wp
    20341928    IF ( dp_external )  THEN
    20351929!
    2036 !--    Set the starting level dp_level_ind_b only if it has not been set before
    2037 !--    (e.g. in init_grid).
     1930!--    Set the starting level dp_level_ind_b only if it has not been set before (e.g. in init_grid).
    20381931       IF ( dp_level_ind_b == 0 )  THEN
    20391932          ind_array = MINLOC( ABS( dp_level_b - zu ) )
    2040           dp_level_ind_b = ind_array(1) - 1 + nzb 
     1933          dp_level_ind_b = ind_array(1) - 1 + nzb
    20411934                                        ! MINLOC uses lower array bound 1
    20421935       ENDIF
     
    20441937          dp_smooth_factor(:dp_level_ind_b) = 0.0_wp
    20451938          DO  k = dp_level_ind_b+1, nzt
    2046              dp_smooth_factor(k) = 0.5_wp * ( 1.0_wp + SIN( pi *               &
    2047                         ( REAL( k - dp_level_ind_b, KIND=wp ) /                &
    2048                           REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5_wp ) ) )
     1939             dp_smooth_factor(k) = 0.5_wp * ( 1.0_wp + SIN( pi *                                   &
     1940                                             ( REAL( k - dp_level_ind_b, KIND=wp ) /               &
     1941                                               REAL( nzt - dp_level_ind_b, KIND=wp ) - 0.5_wp ) ) )
    20491942          ENDDO
    20501943       ENDIF
     
    20521945
    20531946!
    2054 !-- Initialize damping zone for the potential temperature in case of
    2055 !-- non-cyclic lateral boundaries. The damping zone has the maximum value
    2056 !-- at the inflow boundary and decreases to zero at pt_damping_width.
     1947!-- Initialize damping zone for the potential temperature in case of non-cyclic lateral boundaries.
     1948!-- The damping zone has the maximum value at the inflow boundary and decreases to zero at
     1949!-- pt_damping_width.
    20571950    ptdf_x = 0.0_wp
    20581951    ptdf_y = 0.0_wp
     
    20601953       DO  i = nxl, nxr
    20611954          IF ( ( i * dx ) < pt_damping_width )  THEN
    2062              ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp *              &
    2063                             REAL( pt_damping_width - i * dx, KIND=wp ) / (     &
    2064                             REAL( pt_damping_width, KIND=wp ) ) ) )**2
    2065           ENDIF
     1955             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp *                                  &
     1956                                                    REAL( pt_damping_width - i * dx, KIND=wp ) /   &
     1957                                                    REAL( pt_damping_width, KIND=wp ) ) )**2
     1958                                  ENDIF
    20661959       ENDDO
    20671960    ELSEIF ( bc_lr_raddir )  THEN
    20681961       DO  i = nxl, nxr
    20691962          IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) )  THEN
    2070              ptdf_x(i) = pt_damping_factor *                                   &
    2071                          SIN( pi * 0.5_wp *                                    &
    2072                                  ( ( i - nx ) * dx + pt_damping_width ) /      &
    2073                                  REAL( pt_damping_width, KIND=wp ) )**2
    2074           ENDIF
    2075        ENDDO
     1963             ptdf_x(i) = pt_damping_factor * SIN( pi * 0.5_wp *                                    &
     1964                                                  ( ( i - nx ) * dx + pt_damping_width ) /         &
     1965                                                  REAL( pt_damping_width, KIND=wp ) )**2
     1966          ENDIF
     1967       ENDDO
    20761968    ELSEIF ( bc_ns_dirrad )  THEN
    20771969       DO  j = nys, nyn
    20781970          IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) )  THEN
    2079              ptdf_y(j) = pt_damping_factor *                                   &
    2080                          SIN( pi * 0.5_wp *                                    &
    2081                                  ( ( j - ny ) * dy + pt_damping_width ) /      &
    2082                                  REAL( pt_damping_width, KIND=wp ) )**2
    2083           ENDIF
    2084        ENDDO
     1971             ptdf_y(j) = pt_damping_factor * SIN( pi * 0.5_wp *                                    &
     1972                                                  ( ( j - ny ) * dy + pt_damping_width ) /         &
     1973                                                  REAL( pt_damping_width, KIND=wp ) )**2
     1974          ENDIF
     1975       ENDDO
    20851976    ELSEIF ( bc_ns_raddir )  THEN
    20861977       DO  j = nys, nyn
    20871978          IF ( ( j * dy ) < pt_damping_width )  THEN
    2088              ptdf_y(j) = pt_damping_factor *                                   &
    2089                          SIN( pi * 0.5_wp *                                    &
    2090                                 ( pt_damping_width - j * dy ) /                &
    2091                                 REAL( pt_damping_width, KIND=wp ) )**2
     1979             ptdf_y(j) = pt_damping_factor * SIN( pi * 0.5_wp *                                    &
     1980                                                  ( pt_damping_width - j * dy ) /                  &
     1981                                                  REAL( pt_damping_width, KIND=wp ) )**2
    20921982          ENDIF
    20931983       ENDDO
     
    20951985
    20961986!
    2097 !-- Input binary data file is not needed anymore. This line must be placed
    2098 !-- after call of user_init!
     1987!-- Input binary data file is not needed anymore. This line must be placed after call of user_init!
    20991988    CALL close_file( 13 )
    21001989!
    2101 !-- In case of nesting, put an barrier to assure that all parent and child
    2102 !-- domains finished initialization.
     1990!-- In case of nesting, put an barrier to assure that all parent and child domains finished
     1991!-- initialization.
    21031992#if defined( __parallel )
    21041993    IF ( nested_run )  CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
  • palm/trunk/SOURCE/init_advec.f90

    r4360 r4648  
    11!> @file init_advec.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
    98!
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
    1312!
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    2120! -----------------
    22 ! 
    23 ! 
     21!
     22!
    2423! Former revisions:
    2524! -----------------
    2625! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4360 2020-01-07 11:25:50Z suehring
    2729! Corrected "Former revisions" section
    28 ! 
     30!
    2931! 3655 2019-01-07 16:51:22Z knoop
    3032! Corrected "Former revisions" section
     
    3739! ------------
    3840!> Initialize constant coefficients and parameters for certain advection schemes.
    39 !------------------------------------------------------------------------------!
     41!--------------------------------------------------------------------------------------------------!
    4042 SUBROUTINE init_advec
    41  
    4243
    43     USE advection,                                                             &
     44
     45    USE advection,                                                                                 &
    4446        ONLY:  aex, bex, dex, eex
    45        
     47
    4648    USE kinds
    47    
    48     USE control_parameters,                                                    &
     49
     50    USE control_parameters,                                                                        &
    4951        ONLY:  scalar_advec
    5052
     
    5254
    5355    INTEGER(iwp) ::  i          !<
    54     INTEGER(iwp) ::  intervals  !< 
     56    INTEGER(iwp) ::  intervals  !<
    5557    INTEGER(iwp) ::  j          !<
    56    
     58
    5759    REAL(wp) :: delt   !<
    5860    REAL(wp) :: dn     !<
     
    8890             ex1 = dn * EXP( -dn ) - EXP( 0.5_wp * dn ) + EXP( -0.5_wp * dn )
    8991             ex2 = EXP( dn ) - EXP( -dn )
    90              ex3 = EXP( -dn ) * ( 1.0_wp - dn ) - 0.5_wp * EXP(  0.5_wp * dn ) &
     92             ex3 = EXP( -dn ) * ( 1.0_wp - dn ) - 0.5_wp * EXP(  0.5_wp * dn )                     &
    9193                                                - 0.5_wp * EXP( -0.5_wp * dn )
    9294             ex4 = EXP( dn ) + EXP( -dn )
  • palm/trunk/SOURCE/init_coupling.f90

    r4564 r4648  
    11!> @file init_coupling.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
    98!
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
    1312!
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    2120! -----------------
    22 ! 
    23 ! 
     21!
     22!
    2423! Former revisions:
    2524! ------------------
    2625! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4564 2020-06-12 14:03:36Z raasch
    2729! Vertical nesting method of Huq et al. (2019) removed
    28 ! 
     30!
    2931! 4444 2020-03-05 15:59:50Z raasch
    3032! bugfix: cpp-directives for serial mode added
    31 ! 
     33!
    3234! 4360 2020-01-07 11:25:50Z suehring
    3335! Corrected "Former revisions" section
    34 ! 
     36!
    3537! 3655 2019-01-07 16:51:22Z knoop
    3638! references to mrun replaced by palmrun, and updated
     
    4143! Description:
    4244! ------------
    43 !> Initializing coupling via MPI-1 or MPI-2 if the coupled version of PALM is
    44 !> called.
    45 !------------------------------------------------------------------------------!
     45!> Initializing coupling via MPI-1 or MPI-2 if the coupled version of PALM is called.
     46!--------------------------------------------------------------------------------------------------!
    4647  SUBROUTINE init_coupling
    47  
    4848
    49     USE control_parameters,                                                    &
     49
     50    USE control_parameters,                                                                        &
    5051        ONLY:  coupling_char, coupling_mode
    51        
     52
    5253    USE kinds
    53    
     54
    5455    USE pegrid
    5556
     
    6263    INTEGER(iwp) ::  inter_color  !<
    6364#endif
    64    
     65
    6566    INTEGER(iwp), DIMENSION(:) ::  bc_data(0:3) = 0  !<
    6667
    6768!
    68 !-- Get information about the coupling mode from the environment variable
    69 !-- which has been set by the mpiexec command.
    70 !-- This method is currently not used because the mpiexec command is not
    71 !-- available on some machines
     69!-- Get information about the coupling mode from the environment variable which has been set by the
     70!-- mpiexec command.
     71!-- This method is currently not used because the mpiexec command is not available on some machines.
    7272!    CALL GET_ENVIRONMENT_VARIABLE( 'coupling_mode', coupling_mode, i )
    7373!    IF ( i == 0 )  coupling_mode = 'uncoupled'
     
    7575
    7676!
    77 !-- Get information about the coupling mode from standard input (PE0 only) and
    78 !-- distribute it to the other PEs. Distribute PEs to 2 new communicators.
     77!-- Get information about the coupling mode from standard input (PE0 only) and distribute it to the
     78!-- other PEs. Distribute PEs to 2 new communicators.
    7979!-- ATTENTION: numprocs will be reset according to the new communicators
    8080#if defined ( __parallel )
     
    9191
    9292!
    93 !--    Check if '_O' has to be used as file extension in an uncoupled ocean
    94 !--    run. This is required, if this run shall be continued as a coupled run.
     93!--    Check if '_O' has to be used as file extension in an uncoupled ocean run. This is required,
     94!--    if this run shall be continued as a coupled run.
    9595       IF ( TRIM( coupling_mode ) == 'precursor_ocean' )  bc_data(3) = 1
    9696
     
    127127
    128128!
    129 !--    Write a flag file for the ocean model and the other atmosphere
    130 !--    processes.
     129!--    Write a flag file for the ocean model and the other atmosphere processes.
    131130       OPEN( 90, FILE='COUPLING_PORT_OPENED', FORM='FORMATTED' )
    132131       WRITE ( 90, '(''TRUE'')' )
     
    136135
    137136!
    138 !-- In case of a precursor ocean run (followed by a coupled run), or a
    139 !-- coupled atmosphere-ocean run, set the file extension for the ocean files
    140     IF ( TRIM( coupling_mode ) == 'ocean_to_atmosphere' .OR. bc_data(3) == 1 ) &
    141     THEN
     137!-- In case of a precursor ocean run (followed by a coupled run), or a coupled atmosphere-ocean run,
     138!-- set the file extension for the ocean files.
     139    IF ( TRIM( coupling_mode ) == 'ocean_to_atmosphere' .OR. bc_data(3) == 1 )  THEN
    142140       coupling_char = '_O'
    143141    ENDIF
  • palm/trunk/SOURCE/init_grid.f90

    r4630 r4648  
    11!> @file init_grid.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
     
    2524! -----------------
    2625! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4630 2020-07-30 14:54:34Z suehring
    2729! In case of ASCII topography input flag grid points as terrain and building.
    2830!
    2931! 4601 2020-07-14 12:06:09Z suehring
    3032! Minor formatting adjustments
    31 ! 
     33!
    3234! 4564 2020-06-12 14:03:36Z raasch
    3335! Vertical nesting method of Huq et al. (2019) removed
    34 ! 
     36!
    3537! 4543 2020-05-20 14:12:22Z gronemeier
    3638! Remove non-required check for canyon height
    37 ! 
     39!
    3840! 4507 2020-04-22 18:21:45Z gronemeier
    3941! update origin_z with shifting height of orography (oro_min)
    40 ! 
     42!
    4143! 4457 2020-03-11 14:20:43Z raasch
    4244! use statement for exchange horiz added,
    4345! bugfix for call of exchange horiz 2d
    44 ! 
     46!
    4547! 4444 2020-03-05 15:59:50Z raasch
    4648! bugfix: cpp-directives for serial mode added
    47 ! 
     49!
    4850! 4414 2020-02-19 20:16:04Z suehring
    4951! - Remove deprecated topography arrays nzb_s_inner, nzb_u_inner, etc.
    50 ! - Move initialization of boundary conditions and multigrid into an extra
    51 !   module interface.
    52 !
     52! - Move initialization of boundary conditions and multigrid into an extra module interface.
     53!
    5354! 4386 2020-01-27 15:07:30Z Giersch
    54 ! Allocation statements, comments, naming of variables revised and _wp added to
    55 ! real type values
    56 !
     55! Allocation statements, comments, naming of variables revised and _wp added to real type values
     56!
    5757! 4360 2020-01-07 11:25:50Z suehring
    5858! Revise error messages for generic tunnel setup.
    59 ! 
     59!
    6060! 4346 2019-12-18 11:55:56Z motisi
    61 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    62 ! topography information used in wall_flags_static_0
    63 ! 
     61! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     62! information used in wall_flags_static_0
     63!
    6464! 4340 2019-12-16 08:17:03Z Giersch
    6565! Topography closed channel flow with symmetric boundaries implemented
    66 ! 
     66!
    6767! 4329 2019-12-10 15:46:36Z motisi
    6868! Renamed wall_flags_0 to wall_flags_static_0
    69 ! 
     69!
    7070! 4328 2019-12-09 18:53:04Z suehring
    7171! Minor change in nzb_max computation. Commentation added.
    72 ! 
     72!
    7373! 4314 2019-11-29 10:29:20Z suehring
    74 ! Set additional topography flag 4 to mark topography grid points emerged
    75 ! from the filtering process.
    76 !
     74! Set additional topography flag 4 to mark topography grid points emerged from the filtering process.
     75!
    7776! 4294 2019-11-13 18:34:16Z suehring
    78 ! Bugfix, always set bit 5 and 6 of wall_flags, indicating terrain- and
    79 ! building surfaces in all  cases, in order to enable terrain-following output
    80 ! also when no land- or urban-surface model is applied.
    81 ! 
     77! Bugfix, always set bit 5 and 6 of wall_flags, indicating terrain- and building surfaces in all
     78! cases, in order to enable terrain-following output also when no land- or urban-surface model is
     79! applied.
     80!
    8281! 4265 2019-10-15 16:16:24Z suehring
    83 ! Bugfix for last commit, exchange oro_max variable only when it is allocated
    84 ! (not necessarily the case when topography is input from ASCII file).
    85 ! 
     82! Bugfix for last commit, exchange oro_max variable only when it is allocated (not necessarily the
     83! case when topography is input from ASCII file).
     84!
    8685! 4245 2019-09-30 08:40:37Z pavelkrc
    8786! Store oro_max (building z-offset) in 2D for building surfaces
    88 ! 
     87!
    8988! 4189 2019-08-26 16:19:38Z suehring
    9089! - Add check for proper setting of namelist parameter topography
    9190! - Set flag to indicate land surfaces in case no topography is provided
    92 ! 
     91!
    9392! 4182 2019-08-22 15:20:23Z scharf
    9493! Corrected "Former revisions" section
    95 ! 
     94!
    9695! 4168 2019-08-16 13:50:17Z suehring
    97 ! Pre-calculate topography top index and store it on an array (replaces former
    98 ! functions get_topography_top_index)
    99 ! 
     96! Pre-calculate topography top index and store it on an array (replaces former functions
     97! get_topography_top_index)
     98!
    10099! 4159 2019-08-15 13:31:35Z suehring
    101 ! Revision of topography processing. This was not consistent between 2D and 3D
    102 ! buildings.
    103 !
     100! Revision of topography processing. This was not consistent between 2D and 3D buildings.
     101!
    104102! 4144 2019-08-06 09:11:47Z raasch
    105103! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
    106 ! 
     104!
    107105! 4115 2019-07-24 12:50:49Z suehring
    108 ! Bugfix in setting near-surface flag 24, inidicating wall-bounded grid points 
    109 ! 
     106! Bugfix in setting near-surface flag 24, inidicating wall-bounded grid points
     107!
    110108! 4110 2019-07-22 17:05:21Z suehring
    111109! - Separate initialization of advection flags for momentum and scalars.
    112110! - Change subroutine interface for ws_init_flags_scalar to pass boundary flags
    113 ! 
     111!
    114112! 4109 2019-07-22 17:00:34Z suehring
    115113! Fix bad commit
    116 ! 
     114!
    117115! 3926 2019-04-23 12:56:42Z suehring
    118 ! Minor bugfix in building mapping when all building IDs in the model domain
    119 ! are missing
    120 !
     116! Minor bugfix in building mapping when all building IDs in the model domain are missing
     117!
    121118! 3857 2019-04-03 13:00:16Z knoop
    122 ! In projection of non-building 3D objects onto numerical grid remove
    123 ! dependency on building_type
    124 !
     119! In projection of non-building 3D objects onto numerical grid remove dependency on building_type
     120!
    125121! 3763 2019-02-25 17:33:49Z suehring
    126 ! Replace work-around for ghost point exchange of 1-byte arrays with specific
    127 ! routine as already done in other routines
    128 ! 
     122! Replace work-around for ghost point exchange of 1-byte arrays with specific routine as already
     123! done in other routines
     124!
    129125! 3761 2019-02-25 15:31:42Z raasch
    130126! unused variables removed
    131 ! 
     127!
    132128! 3661 2019-01-08 18:22:50Z suehring
    133 ! Remove setting of nzb_max to nzt at non-cyclic boundary PEs, instead,
    134 ! order degradation of advection scheme is handeled directly in advec_ws
    135 ! 
     129! Remove setting of nzb_max to nzt at non-cyclic boundary PEs, instead, order degradation of
     130! advection scheme is handeled directly in advec_ws
     131!
    136132! 3655 2019-01-07 16:51:22Z knoop
    137133! Comment added
     
    142138!
    143139! Description:
    144 ! -----------------------------------------------------------------------------!
     140! -------------------------------------------------------------------------------------------------!
    145141!> Creating grid depending constants
    146142!> @todo: Rearrange topo flag list
    147 !> @todo: reference 3D buildings on top of orography is not tested and may need
    148 !>        further improvement for steep slopes
    149 !> @todo: Use more advanced setting of building type at filled holes 
    150 !------------------------------------------------------------------------------!
     143!> @todo: reference 3D buildings on top of orography is not tested and may need further improvement
     144!>        for steep slopes
     145!> @todo: Use more advanced setting of building type at filled holes
     146!--------------------------------------------------------------------------------------------------!
    151147 SUBROUTINE init_grid
    152148
    153     USE arrays_3d,                                                             &
     149    USE arrays_3d,                                                                                 &
    154150        ONLY:  dd2zu, ddzu, ddzu_pres, ddzw, dzu, dzw, x, xu, y, yv, zu, zw
    155151
    156     USE control_parameters,                                                    &
    157         ONLY:  constant_flux_layer, dz, dz_max, dz_stretch_factor,             &
    158                dz_stretch_factor_array, dz_stretch_level, dz_stretch_level_end,&
    159                dz_stretch_level_end_index, dz_stretch_level_start_index,       &
    160                dz_stretch_level_start, ibc_uv_b, message_string,               &
    161                number_stretch_level_end,                                       &
    162                number_stretch_level_start,                                     &
    163                ocean_mode,                                                     &
    164                psolver,                                                        &
    165                symmetry_flag,                                                  &
    166                topography,                                                     &
     152    USE control_parameters,                                                                        &
     153        ONLY:  constant_flux_layer, dz, dz_max, dz_stretch_factor,                                 &
     154               dz_stretch_factor_array, dz_stretch_level, dz_stretch_level_end,                    &
     155               dz_stretch_level_end_index, dz_stretch_level_start_index,                           &
     156               dz_stretch_level_start, ibc_uv_b, message_string,                                   &
     157               number_stretch_level_end,                                                           &
     158               number_stretch_level_start,                                                         &
     159               ocean_mode,                                                                         &
     160               psolver,                                                                            &
     161               symmetry_flag,                                                                      &
     162               topography,                                                                         &
    167163               use_surface_fluxes
    168164
    169     USE grid_variables,                                                        &
     165    USE grid_variables,                                                                            &
    170166        ONLY:  ddx, ddx2, ddy, ddy2, dx, dx2, dy, dy2, zu_s_inner, zw_w_inner
    171167
    172     USE indices,                                                               &
    173         ONLY:  nbgp,                                                           &
    174                nx,                                                             &
    175                nxl,                                                            &
    176                nxlg,                                                           &
    177                nxr,                                                            &
    178                nxrg,                                                           &
    179                ny,                                                             &
    180                nyn,                                                            &
    181                nyng,                                                           &
    182                nys,                                                            &
    183                nysg,                                                           &
    184                nz,                                                             &
    185                nzb,                                                            &
    186                nzb_diff,                                                       &
    187                nzb_max,                                                        &
    188                nzt,                                                            &
    189                topo_top_ind,                                                   &
     168    USE indices,                                                                                   &
     169        ONLY:  nbgp,                                                                               &
     170               nx,                                                                                 &
     171               nxl,                                                                                &
     172               nxlg,                                                                               &
     173               nxr,                                                                                &
     174               nxrg,                                                                               &
     175               ny,                                                                                 &
     176               nyn,                                                                                &
     177               nyng,                                                                               &
     178               nys,                                                                                &
     179               nysg,                                                                               &
     180               nz,                                                                                 &
     181               nzb,                                                                                &
     182               nzb_diff,                                                                           &
     183               nzb_max,                                                                            &
     184               nzt,                                                                                &
     185               topo_top_ind,                                                                       &
    190186               topo_min_level
    191187
     
    196192    IMPLICIT NONE
    197193
    198     INTEGER(iwp) ::  i             !< index variable along x 
     194    INTEGER(iwp) ::  i             !< index variable along x
    199195    INTEGER(iwp) ::  j             !< index variable along y
    200196    INTEGER(iwp) ::  k             !< index variable along z
    201197    INTEGER(iwp) ::  k_top         !< topography top index on local PE
    202198    INTEGER(iwp) ::  n             !< loop variable for stretching
    203     INTEGER(iwp) ::  number_dz     !< number of user-specified dz values       
     199    INTEGER(iwp) ::  number_dz     !< number of user-specified dz values
    204200    INTEGER(iwp) ::  nzb_local_max !< vertical grid index of maximum topography height
    205201    INTEGER(iwp) ::  nzb_local_min !< vertical grid index of minimum topography height
     
    209205    REAL(wp) ::  dz_level_end  !< distance between calculated height level for u/v-grid and user-specified end level for stretching
    210206    REAL(wp) ::  dz_stretched  !< stretched vertical grid spacing
    211    
    212     REAL(wp), DIMENSION(:), ALLOCATABLE ::  min_dz_stretch_level_end !< Array that contains all minimum heights where the stretching can end
     207
     208    REAL(wp), DIMENSION(:), ALLOCATABLE ::  min_dz_stretch_level_end !< Array that contains all minimum heights where the stretching
     209                                                                     !< can end
    213210
    214211
     
    224221    ALLOCATE( x(0:nx) )
    225222    ALLOCATE( xu(0:nx) )
    226    
     223
    227224    DO i = 0, nx
    228225       xu(i) = i * dx
     
    232229    ALLOCATE( y(0:ny) )
    233230    ALLOCATE( yv(0:ny) )
    234    
     231
    235232    DO j = 0, ny
    236233       yv(j) = j * dy
     
    247244
    248245!
    249 !-- For constructing an appropriate grid, the vertical grid spacing dz has to
    250 !-- be specified with a non-negative value in the parameter file
     246!-- For constructing an appropriate grid, the vertical grid spacing dz has to be specified with a
     247!-- non-negative value in the parameter file.
    251248    IF ( dz(1) == -1.0_wp )  THEN
    252249       message_string = 'missing dz'
    253        CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 ) 
     250       CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 )
    254251    ELSEIF ( dz(1) <= 0.0_wp )  THEN
    255252       WRITE( message_string, * ) 'dz=',dz(1),' <= 0.0'
     
    258255
    259256!
    260 !-- Initialize dz_stretch_level_start with the value of dz_stretch_level
    261 !-- if it was set by the user
     257!-- Initialize dz_stretch_level_start with the value of dz_stretch_level if it was set by the user.
    262258    IF ( dz_stretch_level /= -9999999.9_wp ) THEN
    263259       dz_stretch_level_start(1) = dz_stretch_level
    264260    ENDIF
    265        
    266 !
    267 !-- Determine number of dz values and stretching levels specified by the
    268 !-- user to allow right controlling of the stretching mechanism and to
    269 !-- perform error checks. The additional requirement that dz /= dz_max
    270 !-- for counting number of user-specified dz values is necessary. Otherwise
    271 !-- restarts would abort if the old stretching mechanism with dz_stretch_level
    272 !-- is used (Attention: The user is not allowed to specify a dz value equal
    273 !-- to the default of dz_max = 999.0).
    274     number_dz = COUNT( dz /= -1.0_wp .AND. dz /= dz_max)
    275     number_stretch_level_start = COUNT( dz_stretch_level_start /=              &
    276                                        -9999999.9_wp )
    277     number_stretch_level_end = COUNT( dz_stretch_level_end /=                  &
    278                                       9999999.9_wp )
    279 
    280 !
    281 !-- The number of specified end levels +1 has to be the same as the number
     261
     262!
     263!-- Determine number of dz values and stretching levels specified by the user to allow right
     264!-- controlling of the stretching mechanism and to perform error checks. The additional requirement
     265!-- that dz /= dz_max for counting number of user-specified dz values is necessary. Otherwise
     266!-- restarts would abort if the old stretching mechanism with dz_stretch_level is used (Attention:
     267!-- The user is not allowed to specify a dz value equal to the default of dz_max = 999.0).
     268    number_dz = COUNT( dz /= -1.0_wp  .AND.  dz /= dz_max)
     269    number_stretch_level_start = COUNT( dz_stretch_level_start /= -9999999.9_wp )
     270    number_stretch_level_end = COUNT( dz_stretch_level_end /= 9999999.9_wp )
     271
     272!
     273!-- The number of specified end levels +1 has to be the same as the number
    282274!-- of specified dz values
    283275    IF ( number_dz /= number_stretch_level_end + 1 ) THEN
    284        WRITE( message_string, * ) 'The number of values for dz = ',            &
    285                                    number_dz, 'has to be the same as& ',       &
    286                                    'the number of values for ',                &
    287                                    'dz_stretch_level_end + 1 = ',              &
    288                                    number_stretch_level_end+1
     276       WRITE( message_string, * )  'The number of values for dz = ', number_dz,                    &
     277                                   'has to be the same as& ', 'the number of values for ',         &
     278                                   'dz_stretch_level_end + 1 = ', number_stretch_level_end+1
    289279       CALL message( 'init_grid', 'PA0156', 1, 2, 0, 6, 0 )
    290280    ENDIF
    291    
    292 !
    293 !-- The number of specified start levels has to be the same or one less than
    294 !-- the number of specified dz values
    295     IF ( number_dz /= number_stretch_level_start + 1 .AND.                     &
    296          number_dz /= number_stretch_level_start ) THEN
    297        WRITE( message_string, * ) 'The number of values for dz = ',            &
    298                                    number_dz, 'has to be the same as or one ', &
    299                                    'more than& the number of values for ',     &
    300                                    'dz_stretch_level_start = ',                &
    301                                    number_stretch_level_start
     281
     282!
     283!-- The number of specified start levels has to be the same or one less than the number of specified
     284!-- dz values
     285    IF ( number_dz /= number_stretch_level_start + 1  .AND.                                        &
     286         number_dz /= number_stretch_level_start )  THEN
     287       WRITE( message_string, * )  'The number of values for dz = ', number_dz,                    &
     288                                   'has to be the same as or one ',                                &
     289                                   'more than& the number of values for ',                         &
     290                                   'dz_stretch_level_start = ', number_stretch_level_start
    302291       CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 )
    303292    ENDIF
    304    
    305 !-- The number of specified start levels has to be the same or one more than
    306 !-- the number of specified end levels
    307     IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND.      &
     293
     294!-- The number of specified start levels has to be the same or one more than the number of specified
     295!-- end levels
     296    IF ( number_stretch_level_start /= number_stretch_level_end + 1  .AND.                         &
    308297         number_stretch_level_start /= number_stretch_level_end ) THEN
    309        WRITE( message_string, * ) 'The number of values for ',                 &
    310                                   'dz_stretch_level_start = ',                 &
    311                                    dz_stretch_level_start, 'has to be the ',   &
    312                                    'same or one more than& the number of ',    &
    313                                    'values for dz_stretch_level_end = ',       &
    314                                    number_stretch_level_end
     298       WRITE( message_string, * )  'The number of values for ',                                    &
     299                                   'dz_stretch_level_start = ', dz_stretch_level_start,            &
     300                                   'has to be the ', 'same or one more than& the number of ',      &
     301                                   'values for dz_stretch_level_end = ', number_stretch_level_end
    315302       CALL message( 'init_grid', 'PA0216', 1, 2, 0, 6, 0 )
    316303    ENDIF
     
    318305!
    319306!-- Initialize dz for the free atmosphere with the value of dz_max
    320     IF ( dz(number_stretch_level_start+1) == -1.0_wp .AND.                     &
    321          number_stretch_level_start /= 0 ) THEN
     307    IF ( dz(number_stretch_level_start+1) == -1.0_wp  .AND.  number_stretch_level_start /= 0 ) THEN
    322308       dz(number_stretch_level_start+1) = dz_max
    323309    ENDIF
    324        
    325 !
    326 !-- Initialize the stretching factor if (infinitely) stretching in the free
    327 !-- atmosphere is desired (dz_stretch_level_end was not specified for the
    328 !-- free atmosphere)
    329     IF ( number_stretch_level_start == number_stretch_level_end + 1 ) THEN
    330        dz_stretch_factor_array(number_stretch_level_start) =                   &
    331        dz_stretch_factor
     310
     311!
     312!-- Initialize the stretching factor if (infinitely) stretching in the free atmosphere is desired
     313!-- (dz_stretch_level_end was not specified for the free atmosphere)
     314    IF ( number_stretch_level_start == number_stretch_level_end + 1 )  THEN
     315       dz_stretch_factor_array(number_stretch_level_start) = dz_stretch_factor
    332316    ENDIF
    333    
     317
    334318!
    335319!-- Allocation of arrays for stretching
     
    339323!-- Define the vertical grid levels. Start with atmosphere branch
    340324    IF ( .NOT. ocean_mode )  THEN
    341    
    342 !
    343 !--    The stretching region has to be large enough to allow for a smooth
    344 !--    transition between two different grid spacings. The number 4 is an
    345 !--    empirical value
     325
     326!
     327!--    The stretching region has to be large enough to allow for a smooth transition between two
     328!--    different grid spacings. The number 4 is an empirical value.
    346329       DO n = 1, number_stretch_level_start
    347           min_dz_stretch_level_end(n) = dz_stretch_level_start(n) +            &
    348                                         4 * MAX( dz(n),dz(n+1) )
     330          min_dz_stretch_level_end(n) = dz_stretch_level_start(n) + 4 * MAX( dz(n),dz(n+1) )
    349331       ENDDO
    350332
    351        IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) >      &
    352                  dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
    353           message_string= 'Each dz_stretch_level_end has to be larger ' //     &
    354                           'than its corresponding value for &' //              &
    355                           'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//   &
     333       IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) >                          &
     334                 dz_stretch_level_end(1:number_stretch_level_start) ) )  THEN
     335          message_string= 'Each dz_stretch_level_end has to be larger ' //                         &
     336                          'than its corresponding value for &' //                                  &
     337                          'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//                       &
    356338                          'to allow for smooth grid stretching'
    357339          CALL message( 'init_grid', 'PA0224', 1, 2, 0, 6, 0 )
    358340       ENDIF
    359        
    360 !
    361 !--    Stretching must not be applied within the surface layer
    362 !--    (first two grid points). For the default case dz_stretch_level_start
    363 !--    is negative. Therefore the absolut value is checked here.
     341
     342!
     343!--    Stretching must not be applied within the surface layer (first two grid points). For the
     344!--    default case dz_stretch_level_start is negative. Therefore the absolut value is checked here.
    364345       IF ( ANY( ABS( dz_stretch_level_start ) <= dz(1) * 1.5_wp ) ) THEN
    365           WRITE( message_string, * ) 'Each dz_stretch_level_start has to be ',&
    366                                      'larger than ', dz(1) * 1.5
     346          WRITE( message_string, * )  'Each dz_stretch_level_start has to be ',                    &
     347                                      'larger than ', dz(1) * 1.5
    367348          CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 )
    368349       ENDIF
    369350
    370351!
    371 !--    The stretching has to start and end on a grid level. Therefore
    372 !--    user-specified values are mapped to the next lowest level. The 
    373 !--    calculation of the first level is realized differently just because of
    374 !--    historical reasons (the advanced/new stretching mechanism was realized 
    375 !--    in a way that results don't change if the old parameters
    376 !--    dz_stretch_level, dz_stretch_factor and dz_max are used)
    377        IF ( number_stretch_level_start /= 0 ) THEN
    378           dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) -        &
    379                                             dz(1)/2.0) / dz(1) )               &
     352!--    The stretching has to start and end on a grid level. Therefore user-specified values are
     353!--    mapped to the next lowest level. The calculation of the first level is realized differently
     354!--    just because of historical reasons (the advanced/new stretching mechanism was realized in a
     355!--    way that results don't change if the old parameters dz_stretch_level, dz_stretch_factor and
     356!--    dz_max are used).
     357       IF ( number_stretch_level_start /= 0 )  THEN
     358          dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) - dz(1)/2.0) / dz(1) )       &
    380359                                      * dz(1) + dz(1)/2.0
    381360       ENDIF
    382        
     361
    383362       IF ( number_stretch_level_start > 1 ) THEN
    384363          DO n = 2, number_stretch_level_start
    385              dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) /      &
    386                                               dz(n) ) * dz(n)
    387           ENDDO
    388        ENDIF
    389        
     364             dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) / dz(n) ) * dz(n)
     365          ENDDO
     366       ENDIF
     367
    390368       IF ( number_stretch_level_end /= 0 ) THEN
    391369          DO n = 1, number_stretch_level_end
    392              dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) /          &
    393                                             dz(n+1) ) * dz(n+1)
     370             dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) / dz(n+1) ) * dz(n+1)
    394371          ENDDO
    395372       ENDIF
     
    397374!
    398375!--    Determine stretching factor if necessary
    399        IF ( number_stretch_level_end >= 1 ) THEN 
     376       IF ( number_stretch_level_end >= 1 ) THEN
    400377          CALL calculate_stretching_factor( number_stretch_level_end )
    401378       ENDIF
     
    403380!
    404381!--    Grid for atmosphere with surface at z=0 (k=0, w-grid).
    405 !--    First compute the u- and v-levels. In case of dirichlet bc for u and v
    406 !--    the first u/v- and w-level (k=0) are defined at same height (z=0).
    407 !--    The second u-level (k=1) corresponds to the top of the
    408 !--    surface layer. In case of symmetric boundaries (closed channel flow),
    409 !--    the first grid point is always at z=0.
    410        IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 .OR.                              &
    411             topography == 'closed_channel' ) THEN
     382!--    First compute the u- and v-levels. In case of dirichlet bc for u and v the first u/v- and
     383!--    w-level (k=0) are defined at same height (z=0).
     384!--    The second u-level (k=1) corresponds to the top of the surface layer. In case of symmetric
     385!--    boundaries (closed channel flow), the first grid point is always at z=0.
     386       IF ( ibc_uv_b == 0  .OR.  ibc_uv_b == 2  .OR.  topography == 'closed_channel' )  THEN
    412387          zu(0) = 0.0_wp
    413388       ELSE
    414389          zu(0) = - dz(1) * 0.5_wp
    415390       ENDIF
    416          
     391
    417392       zu(1) =   dz(1) * 0.5_wp
    418        
    419 !
    420 !--    Determine u and v height levels considering the possibility of grid
    421 !--    stretching in several heights.
     393
     394!
     395!--    Determine u and v height levels considering the possibility of grid stretching in several
     396!--    heights.
    422397       n = 1
    423398       dz_stretch_level_start_index = nzt+1
     
    425400       dz_stretched = dz(1)
    426401
    427 !--    The default value of dz_stretch_level_start is negative, thus the first
    428 !--    condition is true even if no stretching shall be applied. Hence, the
    429 !--    second condition is also necessary.
     402!--    The default value of dz_stretch_level_start is negative, thus the first condition is true
     403!--    even if no stretching shall be applied. Hence, the second condition is also necessary.
    430404       DO  k = 2, nzt+1-symmetry_flag
    431           IF ( dz_stretch_level_start(n) <= zu(k-1) .AND.                      &
    432                dz_stretch_level_start(n) /= -9999999.9_wp ) THEN
     405          IF ( dz_stretch_level_start(n) <= zu(k-1)  .AND.                                         &
     406               dz_stretch_level_start(n) /= -9999999.9_wp )  THEN
    433407             dz_stretched = dz_stretched * dz_stretch_factor_array(n)
    434              
    435              IF ( dz(n) > dz(n+1) ) THEN
     408
     409             IF ( dz(n) > dz(n+1) )  THEN
    436410                dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
    437411             ELSE
    438412                dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
    439413             ENDIF
    440              
    441              IF ( dz_stretch_level_start_index(n) == nzt+1 )                         &
    442              dz_stretch_level_start_index(n) = k-1
    443              
    444           ENDIF
    445          
     414
     415             IF ( dz_stretch_level_start_index(n) == nzt+1 )  dz_stretch_level_start_index(n) = k-1
     416
     417          ENDIF
     418
    446419          zu(k) = zu(k-1) + dz_stretched
    447          
    448 !
    449 !--       Make sure that the stretching ends exactly at dz_stretch_level_end 
    450           dz_level_end = ABS( zu(k) - dz_stretch_level_end(n) ) 
    451          
    452           IF ( dz_level_end  < dz(n+1)/3.0 ) THEN
     420
     421!
     422!--       Make sure that the stretching ends exactly at dz_stretch_level_end
     423          dz_level_end = ABS( zu(k) - dz_stretch_level_end(n) )
     424
     425          IF ( dz_level_end  < dz(n+1)/3.0 )  THEN
    453426             zu(k) = dz_stretch_level_end(n)
    454427             dz_stretched = dz(n+1)
    455428             dz_stretch_level_end_index(n) = k
    456              n = n + 1             
     429             n = n + 1
    457430          ENDIF
    458431       ENDDO
    459        
    460 !
    461 !--    If a closed channel flow is simulated, make sure that grid structure is 
    462 !--    the same for both bottom and top boundary. (Hint: Using a different dz
    463 !--    at the bottom and at the top makes no sense due to symmetric boundaries
    464 !--    where dz should be equal. Therefore, different dz at the bottom and top 
    465 !--    causes an abort (see check_parameters).)
    466        IF ( topography == 'closed_channel' ) THEN
     432
     433!
     434!--    If a closed channel flow is simulated, make sure that grid structure is the same for both
     435!--    bottom and top boundary. (Hint: Using a different dz at the bottom and at the top makes no
     436!--    sense due to symmetric boundaries where dz should be equal. Therefore, different dz at the
     437!--    bottom and top causes an abort (see check_parameters).)
     438       IF ( topography == 'closed_channel' )  THEN
    467439          zu(nzt+1) = zu(nzt) + dz(1) * 0.5_wp
    468440       ENDIF
    469441
    470442!
    471 !--    Compute the w-levels. They are always staggered half-way between the
    472 !--    corresponding u-levels. In case of dirichlet bc for u and v at the
    473 !--    ground the first u- and w-level (k=0) are defined at same height (z=0).
    474 !--    Per default, the top w-level is extrapolated linearly. In case of
    475 !--    a closed channel flow, zu(nzt+1) and zw(nzt) must be set explicitely.
    476 !--    (Hint: Using a different dz at the bottom and at the top makes no sense
    477 !--    due to symmetric boundaries where dz should be equal. Therefore,
    478 !--    different dz at the bottom and top causes an abort (see
    479 !--    check_parameters).)
     443!--    Compute the w-levels. They are always staggered half-way between the corresponding u-levels.
     444!--    In case of dirichlet bc for u and v at the ground the first u- and w-level (k=0) are defined
     445!--    at same height (z=0).
     446!--    Per default, the top w-level is extrapolated linearly. In case of a closed channel flow,
     447!--    zu(nzt+1) and zw(nzt) must be set explicitely.
     448!--    (Hint: Using a different dz at the bottom and at the top makes no sense due to symmetric
     449!--    boundaries where dz should be equal. Therefore, different dz at the bottom and top causes an
     450!--    abort (see check_parameters).)
    480451       zw(0) = 0.0_wp
    481452       DO  k = 1, nzt-symmetry_flag
    482453          zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
    483454       ENDDO
    484        IF ( topography == 'closed_channel' ) THEN
     455       IF ( topography == 'closed_channel' )  THEN
    485456          zw(nzt)   = zw(nzt-1) + dz(1)
    486457          zw(nzt+1) = zw(nzt) + dz(1)
     
    492463
    493464!
    494 !--    The stretching region has to be large enough to allow for a smooth
    495 !--    transition between two different grid spacings. The number 4 is an
    496 !--    empirical value
     465!--    The stretching region has to be large enough to allow for a smooth transition between two
     466!--    different grid spacings. The number 4 is an empirical value
    497467       DO n = 1, number_stretch_level_start
    498           min_dz_stretch_level_end(n) = dz_stretch_level_start(n) -            &
    499                                         4 * MAX( dz(n),dz(n+1) )
     468          min_dz_stretch_level_end(n) = dz_stretch_level_start(n) - 4 * MAX( dz(n),dz(n+1) )
    500469       ENDDO
    501        
    502        IF ( ANY( min_dz_stretch_level_end (1:number_stretch_level_start) <     &
    503                  dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
    504              message_string= 'Each dz_stretch_level_end has to be less ' //   &
    505                              'than its corresponding value for &' //           &
    506                              'dz_stretch_level_start - 4*MAX(dz(n),dz(n+1)) '//&
     470
     471       IF ( ANY( min_dz_stretch_level_end (1:number_stretch_level_start) <                         &
     472                 dz_stretch_level_end(1:number_stretch_level_start) ) )  THEN
     473             message_string= 'Each dz_stretch_level_end has to be less ' //                        &
     474                             'than its corresponding value for &' //                               &
     475                             'dz_stretch_level_start - 4*MAX(dz(n),dz(n+1)) '//                    &
    507476                             'to allow for smooth grid stretching'
    508477             CALL message( 'init_grid', 'PA0224', 1, 2, 0, 6, 0 )
    509478       ENDIF
    510        
    511 !
    512 !--    Stretching must not be applied close to the surface (last two grid
    513 !--    points). For the default case dz_stretch_level_start is negative.
    514        IF ( ANY( dz_stretch_level_start >= - dz(1) * 1.5_wp ) ) THEN
    515           WRITE( message_string, * ) 'Each dz_stretch_level_start has to be ',&
    516                                      'less than ', -dz(1) * 1.5
     479
     480!
     481!--    Stretching must not be applied close to the surface (last two grid points). For the default
     482!--    case dz_stretch_level_start is negative.
     483       IF ( ANY( dz_stretch_level_start >= - dz(1) * 1.5_wp ) )  THEN
     484          WRITE( message_string, * )  'Each dz_stretch_level_start has to be ',                    &
     485                                      'less than ', -dz(1) * 1.5
    517486             CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 )
    518487       ENDIF
    519488
    520489!
    521 !--    The stretching has to start and end on a grid level. Therefore
    522 !--    user-specified values are mapped to the next highest level. The 
    523 !--    calculation of the first level is realized differently just because of
    524 !--    historical reasons (the advanced/new stretching mechanism was realized 
    525 !--    in a way that results don't change if the old parameters
    526 !--    dz_stretch_level, dz_stretch_factor and dz_max are used)
    527        IF ( number_stretch_level_start /= 0 ) THEN
    528           dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) +        &
    529                                             dz(1)/2.0) / dz(1) )               &
     490!--    The stretching has to start and end on a grid level. Therefore user-specified values are
     491!--    mapped to the next highest level. The calculation of the first level is realized differently
     492!--    just because of historical reasons (the advanced/new stretching mechanism was realized in a
     493!--    way that results don't change if the old parameters dz_stretch_level, dz_stretch_factor and
     494!--    dz_max are used)
     495       IF ( number_stretch_level_start /= 0 )  THEN
     496          dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) + dz(1)/2.0) / dz(1) )       &
    530497                                      * dz(1) - dz(1)/2.0
    531498       ENDIF
    532        
    533        IF ( number_stretch_level_start > 1 ) THEN
     499
     500       IF ( number_stretch_level_start > 1 )  THEN
    534501          DO n = 2, number_stretch_level_start
    535              dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) /      &
    536                                               dz(n) ) * dz(n)
    537           ENDDO
    538        ENDIF
    539        
    540        IF ( number_stretch_level_end /= 0 ) THEN
     502             dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) / dz(n) ) * dz(n)
     503          ENDDO
     504       ENDIF
     505
     506       IF ( number_stretch_level_end /= 0 )  THEN
    541507          DO n = 1, number_stretch_level_end
    542              dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) /          &
    543                                             dz(n+1) ) * dz(n+1)
    544           ENDDO
    545        ENDIF
    546        
     508             dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) / dz(n+1) ) * dz(n+1)
     509          ENDDO
     510       ENDIF
     511
    547512!
    548513!--    Determine stretching factor if necessary
    549        IF ( number_stretch_level_end >= 1 ) THEN 
     514       IF ( number_stretch_level_end >= 1 ) THEN
    550515          CALL calculate_stretching_factor( number_stretch_level_end )
    551516       ENDIF
     
    553518!
    554519!--    Grid for ocean with free water surface is at k=nzt (w-grid).
    555 !--    In case of neumann bc at the ground the first first u-level (k=0) lies
    556 !--    below the first w-level (k=0). In case of dirichlet bc the first u- and
    557 !--    w-level are defined at same height, but staggered from the second level.
     520!--    In case of neumann bc at the ground the first first u-level (k=0) lies below the first
     521!--    w-level (k=0). In case of dirichlet bc the first u- and w-level are defined at same height,
     522!--    but staggered from the second level.
    558523!--    The second u-level (k=1) corresponds to the top of the surface layer.
    559524!--    z values are negative starting from z=0 (surface)
     
    562527
    563528!
    564 !--    Determine u and v height levels considering the possibility of grid
    565 !--    stretching in several heights.
     529!--    Determine u and v height levels considering the possibility of grid stretching in several
     530!--    heights.
    566531       n = 1
    567532       dz_stretch_level_start_index = 0
     
    570535
    571536       DO  k = nzt-1, 0, -1
    572          
    573           IF ( dz_stretch_level_start(n) >= zu(k+1) ) THEN
     537
     538          IF ( dz_stretch_level_start(n) >= zu(k+1) )  THEN
    574539             dz_stretched = dz_stretched * dz_stretch_factor_array(n)
    575540
    576              IF ( dz(n) > dz(n+1) ) THEN
     541             IF ( dz(n) > dz(n+1) )  THEN
    577542                dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
    578543             ELSE
    579544                dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
    580545             ENDIF
    581              
    582              IF ( dz_stretch_level_start_index(n) == 0 )                             &
    583              dz_stretch_level_start_index(n) = k+1
    584              
    585           ENDIF
    586          
     546
     547             IF ( dz_stretch_level_start_index(n) == 0 )  dz_stretch_level_start_index(n) = k+1
     548
     549          ENDIF
     550
    587551          zu(k) = zu(k+1) - dz_stretched
    588          
    589 !
    590 !--       Make sure that the stretching ends exactly at dz_stretch_level_end 
    591           dz_level_end = ABS( zu(k) - dz_stretch_level_end(n) ) 
    592          
    593           IF ( dz_level_end  < dz(n+1)/3.0 ) THEN
     552
     553!
     554!--       Make sure that the stretching ends exactly at dz_stretch_level_end
     555          dz_level_end = ABS( zu(k) - dz_stretch_level_end(n) )
     556
     557          IF ( dz_level_end  < dz(n+1)/3.0 )  THEN
    594558             zu(k) = dz_stretch_level_end(n)
    595559             dz_stretched = dz(n+1)
    596560             dz_stretch_level_end_index(n) = k
    597              n = n + 1             
     561             n = n + 1
    598562          ENDIF
    599563       ENDDO
    600        
    601 !
    602 !--    Compute the w-levels. They are always staggered half-way between the
    603 !--    corresponding u-levels, except in case of dirichlet bc for u and v
    604 !--    at the ground. In this case the first u- and w-level are defined at
    605 !--    same height. The top w-level (nzt+1) is not used but set for
     564
     565!
     566!--    Compute the w-levels. They are always staggered half-way between the corresponding u-levels,
     567!--    except in case of dirichlet bc for u and v at the ground. In this case the first u- and
     568!--    w-level are defined at same height. The top w-level (nzt+1) is not used but set for
    606569!--    consistency, since w and all scalar variables are defined up tp nzt+1.
    607570       zw(nzt+1) = dz(1)
     
    612575
    613576!
    614 !--    In case of dirichlet bc for u and v the first u- and w-level are defined
    615 !--    at same height.
     577!--    In case of dirichlet bc for u and v the first u- and w-level are defined at same height.
    616578       IF ( ibc_uv_b == 0 ) THEN
    617579          zu(0) = zw(0)
     
    632594       dd2zu(k) = 1.0_wp / ( dzu(k) + dzu(k+1) )
    633595    ENDDO
    634    
    635 !   
    636 !-- The FFT- SOR-pressure solvers assume grid spacings of a staggered grid
    637 !-- everywhere. For the actual grid, the grid spacing at the lowest level
    638 !-- is only dz/2, but should be dz. Therefore, an additional array
    639 !-- containing with appropriate grid information is created for these
    640 !-- solvers.
     596
     597!
     598!-- The FFT- SOR-pressure solvers assume grid spacings of a staggered grid everywhere. For the
     599!-- actual grid, the grid spacing at the lowest level is only dz/2, but should be dz. Therefore, an
     600!-- additional array containing with appropriate grid information is created for these solvers.
    641601    IF ( psolver(1:9) /= 'multigrid' )  THEN
    642602       ALLOCATE( ddzu_pres(1:nzt+1) )
     
    659619    topo = 0
    660620!
    661 !-- Initialize topography by generic topography or read topography from file. 
     621!-- Initialize topography by generic topography or read topography from file.
    662622    CALL init_topo( topo )
    663623!
    664 !-- Set flags to mask topography on the grid. 
     624!-- Set flags to mask topography on the grid.
    665625    CALL set_topo_flags( topo )
    666626
    667627!
    668 !-- Determine the maximum level of topography. It is used for
    669 !-- steering the degradation of order of the applied advection scheme,
    670 !-- as well in the lpm.
     628!-- Determine the maximum level of topography. It is used for steering the degradation of order of
     629!-- the applied advection scheme, as well in the lpm.
    671630    k_top = 0
    672631    DO  i = nxl, nxr
     
    678637    ENDDO
    679638#if defined( __parallel )
    680     CALL MPI_ALLREDUCE( k_top, nzb_max, 1, MPI_INTEGER,                        &
    681                         MPI_MAX, comm2d, ierr )
     639    CALL MPI_ALLREDUCE( k_top, nzb_max, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
    682640#else
    683641    nzb_max = k_top
     
    686644!-- Increment nzb_max by 1 in order to allow for proper diverengence correction.
    687645!-- Further, in case topography extents up to the model top, limit to nzt.
    688     nzb_max = MIN( nzb_max+1, nzt ) 
    689 !
    690 !-- Determine minimum index of topography. Usually, this will be nzb. In case
    691 !-- there is elevated topography, however, the lowest topography will be higher.
    692 !-- This index is e.g. used to calculate mean first-grid point atmosphere
    693 !-- temperature, surface pressure and density, etc. .
     646    nzb_max = MIN( nzb_max+1, nzt )
     647!
     648!-- Determine minimum index of topography. Usually, this will be nzb. In case there is elevated
     649!-- topography, however, the lowest topography will be higher.
     650!-- This index is e.g. used to calculate mean first-grid point atmosphere temperature, surface
     651!-- pressure and density, etc. .
    694652    topo_min_level   = 0
    695653#if defined( __parallel )
    696     CALL MPI_ALLREDUCE( MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ),             &
    697                         topo_min_level, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )
     654    CALL MPI_ALLREDUCE( MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ), topo_min_level, 1, MPI_INTEGER, &
     655                        MPI_MIN, comm2d, ierr )
    698656#else
    699657    topo_min_level = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
     
    701659
    702660!
    703 !-- Check topography for consistency with model domain. Therefore, use
    704 !-- maximum and minium topography-top indices. Note, minimum topography top
    705 !-- index is already calculated. 
     661!-- Check topography for consistency with model domain. Therefore, use maximum and minium
     662!-- topography-top indices. Note, minimum topography top index is already calculated.
    706663    IF ( TRIM( topography ) /= 'flat' )  THEN
    707664#if defined( __parallel )
    708        CALL MPI_ALLREDUCE( MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ),          &
    709                            nzb_local_max, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )               
     665       CALL MPI_ALLREDUCE( MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ), nzb_local_max, 1,            &
     666                           MPI_INTEGER, MPI_MAX, comm2d, ierr )
    710667#else
    711668       nzb_local_max = MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
     
    715672!--    Consistency checks
    716673       IF ( nzb_local_min < 0  .OR.  nzb_local_max  > nz + 1 )  THEN
    717           WRITE( message_string, * ) 'nzb_local values are outside the',       &
    718                                 ' model domain',                               &
    719                                 '&MINVAL( nzb_local ) = ', nzb_local_min,      &
    720                                 '&MAXVAL( nzb_local ) = ', nzb_local_max
     674          WRITE( message_string, * ) 'nzb_local values are outside the model domain',              &
     675                                     '&MINVAL( nzb_local ) = ', nzb_local_min,                     &
     676                                     '&MAXVAL( nzb_local ) = ', nzb_local_max
    721677          CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 )
    722678       ENDIF
    723679    ENDIF
    724680!
    725 !-- Define vertical gridpoint from (or to) which on the usual finite difference
    726 !-- form (which does not use surface fluxes) is applied
     681!-- Define vertical gridpoint from (or to) which on the usual finite difference form (which does not
     682!-- use surface fluxes) is applied.
    727683    IF ( constant_flux_layer  .OR.  use_surface_fluxes )  THEN
    728684       nzb_diff = nzb + 2
     
    733689    IF ( TRIM( topography ) /= 'flat' )  THEN
    734690!
    735 !--    Allocate and set the arrays containing the topography height (for output
    736 !--    reasons only).
     691!--    Allocate and set the arrays containing the topography height (for output reasons only).
    737692       IF ( nxr == nx  .AND.  nyn /= ny )  THEN
    738           ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn),                             &
    739                     zw_w_inner(nxl:nxr+1,nys:nyn) )
     693          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn), zw_w_inner(nxl:nxr+1,nys:nyn) )
    740694       ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
    741           ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1),                             &
    742                     zw_w_inner(nxl:nxr,nys:nyn+1) )
     695          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1), zw_w_inner(nxl:nxr,nys:nyn+1) )
    743696       ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
    744           ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1),                           &
    745                     zw_w_inner(nxl:nxr+1,nys:nyn+1) )
     697          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1), zw_w_inner(nxl:nxr+1,nys:nyn+1) )
    746698       ELSE
    747           ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn),                               &
    748                     zw_w_inner(nxl:nxr,nys:nyn) )
     699          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn), zw_w_inner(nxl:nxr,nys:nyn) )
    749700       ENDIF
    750701
     
    752703       zw_w_inner   = 0.0_wp
    753704!
    754 !--    Determine local topography height on scalar and w-grid. Note, setting
    755 !--    lateral boundary values is not necessary, realized via wall_flags_static_0
    756 !--    array. Further, please note that loop bounds are different from
    757 !--    nxl to nxr and nys to nyn on south and right model boundary, hence,
     705!--    Determine local topography height on scalar and w-grid. Note, setting lateral boundary values
     706!--    is not necessary, realized via wall_flags_static_0 array. Further, please note that loop
     707!--    bounds are different from nxl to nxr and nys to nyn on south and right model boundary, hence,
    758708!--    use intrinsic lbound and ubound functions to infer array bounds.
    759709       DO  i = LBOUND(zu_s_inner, 1), UBOUND(zu_s_inner, 1)
    760710          DO  j = LBOUND(zu_s_inner, 2), UBOUND(zu_s_inner, 2)
    761711!
    762 !--          Topography height on scalar grid. Therefore, determine index of
    763 !--          upward-facing surface element on scalar grid.
     712!--          Topography height on scalar grid. Therefore, determine index of upward-facing surface
     713!--          element on scalar grid.
    764714             zu_s_inner(i,j) = zu(topo_top_ind(j,i,0))
    765715!
    766 !--          Topography height on w grid. Therefore, determine index of
    767 !--          upward-facing surface element on w grid.
     716!--          Topography height on w grid. Therefore, determine index of upward-facing surface
     717!--          element on w grid.
    768718             zw_w_inner(i,j) = zw(topo_top_ind(j,i,3))
    769719          ENDDO
     
    775725
    776726! Description:
    777 ! -----------------------------------------------------------------------------!
    778 !> Calculation of the stretching factor through an iterative method. Ideas were
    779 !> taken from the paper "Regional stretched grid generation and its application
    780 !> to the NCAR RegCM (1999)". Normally, no analytic solution exists because the
    781 !> system of equations has two variables (r,l) but four requirements
    782 !> (l=integer, r=[0,88;1,2], Eq(6), Eq(5) starting from index j=1) which
    783 !> results into an overdetermined system.
    784 !------------------------------------------------------------------------------!
     727! -------------------------------------------------------------------------------------------------!
     728!> Calculation of the stretching factor through an iterative method. Ideas were taken from the paper
     729!> "Regional stretched grid generation and its application to the NCAR RegCM (1999)". Normally, no
     730!> analytic solution exists because the system of equations has two variables (r,l) but four
     731!> requirements  (l=integer, r=[0,88;1,2], Eq(6), Eq(5) starting from index j=1) which results into
     732!> an overdetermined system.
     733!--------------------------------------------------------------------------------------------------!
    785734 SUBROUTINE calculate_stretching_factor( number_end )
    786  
    787     USE control_parameters,                                                    &
    788         ONLY:  dz, dz_stretch_factor_array,                 &
    789                dz_stretch_level_end, dz_stretch_level_start, message_string
    790  
     735
     736    USE control_parameters,                                                                        &
     737        ONLY:  dz, dz_stretch_factor_array, dz_stretch_level_end, dz_stretch_level_start,          &
     738               message_string
     739
    791740    USE kinds
    792    
     741
    793742    IMPLICIT NONE
    794    
    795     INTEGER(iwp) ::  iterations  !< number of iterations until stretch_factor_lower/upper_limit is reached 
    796     INTEGER(iwp) ::  l_rounded   !< after l_rounded grid levels dz(n) is strechted to dz(n+1) with stretch_factor_2
     743
     744    REAL(wp), PARAMETER ::  stretch_factor_interval = 1.0E-06_wp  !< interval for sampling possible stretching factors
     745    REAL(wp), PARAMETER ::  stretch_factor_lower_limit = 0.88_wp  !< lowest possible stretching factor
     746    REAL(wp), PARAMETER ::  stretch_factor_upper_limit = 1.12_wp  !< highest possible stretching factor
     747
     748    INTEGER(iwp) ::  iterations  !< number of iterations until stretch_factor_lower/upper_limit is reached
     749    INTEGER(iwp) ::  l_rounded   !< after l_rounded grid levels dz(n) is strechted to dz(n+1) with stretch_factor_2
    797750    INTEGER(iwp) ::  n           !< loop variable for stretching
    798    
     751
    799752    INTEGER(iwp), INTENT(IN) ::  number_end !< number of user-specified end levels for stretching
    800        
     753
    801754    REAL(wp) ::  delta_l               !< absolute difference between l and l_rounded
    802755    REAL(wp) ::  delta_stretch_factor  !< absolute difference between stretch_factor_1 and stretch_factor_2
    803     REAL(wp) ::  delta_total_new       !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible)
    804     REAL(wp) ::  delta_total_old       !< sum of delta_l and delta_stretch_factor for the last iteration
     756    REAL(wp) ::  delta_total_new       !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as
     757                                       !< possible)
     758    REAL(wp) ::  delta_total_old       !< sum of delta_l and delta_stretch_factor for the last iteration
    805759    REAL(wp) ::  distance              !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region)
    806     REAL(wp) ::  l                     !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly
     760    REAL(wp) ::  l                     !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1
     761                                       !< exactly
    807762    REAL(wp) ::  numerator             !< numerator of the quotient
    808763    REAL(wp) ::  stretch_factor_1      !< stretching factor that fulfil Eq. (5) togehter with l exactly
    809764    REAL(wp) ::  stretch_factor_2      !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly
    810    
    811     REAL(wp) ::  dz_stretch_factor_array_2(9) = 1.08_wp  !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1
    812    
    813     REAL(wp), PARAMETER ::  stretch_factor_interval = 1.0E-06_wp  !< interval for sampling possible stretching factors
    814     REAL(wp), PARAMETER ::  stretch_factor_lower_limit = 0.88_wp  !< lowest possible stretching factor
    815     REAL(wp), PARAMETER ::  stretch_factor_upper_limit = 1.12_wp  !< highest possible stretching factor
    816  
    817  
     765
     766    REAL(wp) ::  dz_stretch_factor_array_2(9) = 1.08_wp  !< Array that contains all stretch_factor_2 that belongs to
     767                                                         !< stretch_factor_1
     768
     769
    818770    l = 0
    819771    DO  n = 1, number_end
    820    
     772
    821773       iterations = 1
    822        stretch_factor_1 = 1.0_wp 
     774       stretch_factor_1 = 1.0_wp
    823775       stretch_factor_2 = 1.0_wp
    824776       delta_total_old = 1.0_wp
    825        
     777
    826778!
    827779!--    First branch for stretching from rough to fine
    828        IF ( dz(n) > dz(n+1) ) THEN
    829           DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit ) 
    830              
     780       IF ( dz(n) > dz(n+1) )  THEN
     781          DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit )
     782
    831783             stretch_factor_1 = 1.0_wp - iterations * stretch_factor_interval
    832              distance = ABS( dz_stretch_level_end(n) -                         &
    833                         dz_stretch_level_start(n) )   
    834              numerator = distance*stretch_factor_1/dz(n) +                     &
    835                          stretch_factor_1 - distance/dz(n)
    836              
    837              IF ( numerator > 0.0_wp ) THEN
     784             distance = ABS( dz_stretch_level_end(n) - dz_stretch_level_start(n) )
     785             numerator = distance * stretch_factor_1 / dz(n) + stretch_factor_1 - distance / dz(n)
     786
     787             IF ( numerator > 0.0_wp )  THEN
    838788                l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0_wp
    839789                l_rounded = NINT( l )
    840790                delta_l = ABS( l_rounded - l ) / l
    841791             ENDIF
    842              
     792
    843793             stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
    844              
    845              delta_stretch_factor = ABS( stretch_factor_1 -                    &
    846                                          stretch_factor_2 ) /                  &
    847                                     stretch_factor_2
    848              
     794
     795             delta_stretch_factor = ABS( stretch_factor_1 - stretch_factor_2 ) / stretch_factor_2
     796
    849797             delta_total_new = delta_l + delta_stretch_factor
    850798
    851799!
    852 !--          stretch_factor_1 is taken to guarantee that the stretching
    853 !--          procedure ends as close as possible to dz_stretch_level_end.
    854 !--          stretch_factor_2 would guarantee that the stretched dz(n) is
    855 !--          equal to dz(n+1) after l_rounded grid levels.
    856              IF (delta_total_new < delta_total_old) THEN
     800!--          stretch_factor_1 is taken to guarantee that the stretching procedure ends as close as
     801!--          possible to dz_stretch_level_end.
     802!--          stretch_factor_2 would guarantee that the stretched dz(n) is equal to dz(n+1) after
     803!--          l_rounded grid levels.
     804             IF (delta_total_new < delta_total_old)  THEN
    857805                dz_stretch_factor_array(n) = stretch_factor_1
    858806                dz_stretch_factor_array_2(n) = stretch_factor_2
    859807                delta_total_old = delta_total_new
    860808             ENDIF
    861              
     809
    862810             iterations = iterations + 1
    863            
     811
    864812          ENDDO
    865813
    866814!
    867815!--    Second branch for stretching from fine to rough
    868        ELSEIF ( dz(n) < dz(n+1) ) THEN
     816       ELSEIF ( dz(n) < dz(n+1) )  THEN
    869817          DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit )
    870                      
     818
    871819             stretch_factor_1 = 1.0_wp + iterations * stretch_factor_interval
    872              distance = ABS( dz_stretch_level_end(n) -                         &
    873                         dz_stretch_level_start(n) )
    874              numerator = distance*stretch_factor_1/dz(n) +                     &
    875                          stretch_factor_1 - distance/dz(n)
    876              
     820             distance = ABS( dz_stretch_level_end(n) - dz_stretch_level_start(n) )
     821             numerator = distance * stretch_factor_1 / dz(n) + stretch_factor_1 - distance / dz(n)
     822
    877823             l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0_wp
    878824             l_rounded = NINT( l )
    879825             delta_l = ABS( l_rounded - l ) / l
    880              
     826
    881827             stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
    882828
    883              delta_stretch_factor = ABS( stretch_factor_1 -                    &
    884                                         stretch_factor_2 ) /                   &
    885                                         stretch_factor_2
    886              
     829             delta_stretch_factor = ABS( stretch_factor_1 - stretch_factor_2 ) / stretch_factor_2
     830
    887831             delta_total_new = delta_l + delta_stretch_factor
    888              
    889 !
    890 !--          stretch_factor_1 is taken to guarantee that the stretching
    891 !--          procedure ends as close as possible to dz_stretch_level_end.
    892 !--          stretch_factor_2 would guarantee that the stretched dz(n) is
    893 !--          equal to dz(n+1) after l_rounded grid levels.
    894              IF (delta_total_new < delta_total_old) THEN
     832
     833!
     834!--          stretch_factor_1 is taken to guarantee that the stretching procedure ends as close as
     835!--          possible to dz_stretch_level_end.
     836!--          stretch_factor_2 would guarantee that the stretched dz(n) is equal to dz(n+1) after
     837!--          l_rounded grid levels.
     838             IF (delta_total_new < delta_total_old)  THEN
    895839                dz_stretch_factor_array(n) = stretch_factor_1
    896840                dz_stretch_factor_array_2(n) = stretch_factor_2
    897841                delta_total_old = delta_total_new
    898842             ENDIF
    899              
     843
    900844             iterations = iterations + 1
    901845          ENDDO
    902          
     846
    903847       ELSE
    904848          message_string= 'Two adjacent values of dz must be different'
    905849          CALL message( 'init_grid', 'PA0228', 1, 2, 0, 6, 0 )
    906          
    907        ENDIF
    908 
    909 !
    910 !--    Check if also the second stretching factor fits into the allowed
    911 !--    interval. If not, print a warning for the user.
    912        IF ( dz_stretch_factor_array_2(n) < stretch_factor_lower_limit .OR.     &
    913             dz_stretch_factor_array_2(n) > stretch_factor_upper_limit ) THEN
    914           WRITE( message_string, * ) 'stretch_factor_2 = ',                    &
    915                                      dz_stretch_factor_array_2(n), ' which is',&
    916                                      ' responsible for exactly reaching& dz =',&
    917                                       dz(n+1), 'after a specific amount of',   &
    918                                      ' grid levels& exceeds the upper',        &
    919                                      ' limit =', stretch_factor_upper_limit,   &
    920                                      ' &or lower limit = ',                    &
    921                                      stretch_factor_lower_limit
     850
     851       ENDIF
     852
     853!
     854!--    Check if also the second stretching factor fits into the allowed interval. If not, print a
     855!--    warning for the user.
     856       IF ( dz_stretch_factor_array_2(n) < stretch_factor_lower_limit  .OR.                        &
     857            dz_stretch_factor_array_2(n) > stretch_factor_upper_limit )  THEN
     858          WRITE( message_string, * ) 'stretch_factor_2 = ', dz_stretch_factor_array_2(n),          &
     859                                     ' which is', ' responsible for exactly reaching& dz =',       &
     860                                      dz(n+1), 'after a specific amount of',                       &
     861                                     ' grid levels& exceeds the upper',                            &
     862                                     ' limit =', stretch_factor_upper_limit,                       &
     863                                     ' &or lower limit = ', stretch_factor_lower_limit
    922864          CALL message( 'init_grid', 'PA0499', 0, 1, 0, 6, 0 )
    923            
     865
    924866       ENDIF
    925867    ENDDO
    926        
     868
    927869 END SUBROUTINE calculate_stretching_factor
    928  
    929  
     870
     871
    930872! Description:
    931 ! -----------------------------------------------------------------------------!
    932 !> Set temporary topography flags and reference buildings on top of underlying
    933 !> orography.
    934 !------------------------------------------------------------------------------!
     873! -------------------------------------------------------------------------------------------------!
     874!> Set temporary topography flags and reference buildings on top of underlying orography.
     875!--------------------------------------------------------------------------------------------------!
    935876 SUBROUTINE process_topography( topo_3d )
    936877
    937     USE arrays_3d,                                                             &
     878    USE arrays_3d,                                                                                 &
    938879        ONLY:  zu, zw
    939880
    940     USE control_parameters,                                                    &
     881    USE control_parameters,                                                                        &
    941882        ONLY:  bc_lr_cyc, bc_ns_cyc, ocean_mode
    942883
    943     USE exchange_horiz_mod,                                                    &
    944         ONLY:  exchange_horiz_int, exchange_horiz_2d
    945 
    946     USE indices,                                                               &
    947         ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
    948                nzt
    949 
    950     USE netcdf_data_input_mod,                                                 &
    951         ONLY:  buildings_f, building_id_f, building_type_f,                    &
    952                init_model,                                                     &
    953                input_pids_static,                                              &
     884    USE exchange_horiz_mod,                                                                        &
     885        ONLY:  exchange_horiz_2d, exchange_horiz_int
     886
     887    USE indices,                                                                                   &
     888        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nzt
     889
     890    USE netcdf_data_input_mod,                                                                     &
     891        ONLY:  buildings_f, building_id_f, building_type_f,                                        &
     892               init_model,                                                                         &
     893               input_pids_static,                                                                  &
    954894               terrain_height_f
    955895
     
    972912#endif
    973913    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids           !< building IDs on entire model domain
    974     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_final     !< building IDs on entire model domain, multiple occurences are sorted out
     914    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_final     !< building IDs on entire model domain, multiple occurences are
     915                                                                    !< sorted out
    975916    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_final_tmp !< temporary array used for resizing
    976917    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_l         !< building IDs on local subdomain
     
    980921    INTEGER(iwp), DIMENSION(0:numprocs-1) ::  num_buildings_l   !< number of buildings with different ID on local subdomain
    981922
    982     INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo_3d !< input array for 3D topography and dummy array for setting "outer"-flags
    983 
    984     REAL(wp)                            ::  ocean_offset        !< offset to consider inverse vertical coordinate at topography definition
    985     REAL(wp)                            ::  oro_min = 0.0_wp    !< minimum terrain height in entire model domain, used to reference terrain to zero
     923    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo_3d !< input array for 3D topography and dummy array for setting
     924                                                                       !< "outer"-flags
     925
     926    REAL(wp)                            ::  ocean_offset        !< offset to consider inverse vertical coordinate at topography
     927                                                                !< definition
     928    REAL(wp)                            ::  oro_min = 0.0_wp    !< minimum terrain height in entire model domain, used to reference
     929                                                                !< terrain to zero
    986930    REAL(wp), DIMENSION(:), ALLOCATABLE ::  oro_max             !< maximum terrain height occupied by an building with certain id
    987     REAL(wp), DIMENSION(:), ALLOCATABLE ::  oro_max_l           !< maximum terrain height occupied by an building with certain id, on local subdomain
    988 
    989 !
    990 !-- Reference lowest terrain height to zero. This ensures that first,
    991 !-- non-required gird levels (those which lie entirely below the minimum
    992 !-- orography) are avoided, and second, that also negative orography can be used
    993 !-- within the input file.
    994 !-- Please note, in case of a nested run, the global minimum from all parent and
    995 !-- childs need to be remove to avoid steep edges at the child-domain boundaries.
     931    REAL(wp), DIMENSION(:), ALLOCATABLE ::  oro_max_l           !< maximum terrain height occupied by an building with certain id,
     932                                                                !< on local subdomain
     933
     934!
     935!-- Reference lowest terrain height to zero. This ensures that first, non-required gird levels
     936!-- (those which lie entirely below the minimum orography) are avoided, and second, that also
     937!-- negative orography can be used within the input file.
     938!-- Please note, in case of a nested run, the global minimum from all parent and childs needs to be
     939!-- removed to avoid steep edges at the child-domain boundaries.
    996940    IF ( input_pids_static )  THEN
    997    
    998 #if defined( __parallel ) 
    999        CALL MPI_ALLREDUCE( MINVAL( terrain_height_f%var ), oro_min, 1,         &
    1000                            MPI_REAL, MPI_MIN, MPI_COMM_WORLD, ierr )
     941
     942#if defined( __parallel )
     943       CALL MPI_ALLREDUCE( MINVAL( terrain_height_f%var ), oro_min, 1, MPI_REAL, MPI_MIN,          &
     944                           MPI_COMM_WORLD, ierr )
    1001945#else
    1002946       oro_min = MINVAL( terrain_height_f%var )
     
    1007951       init_model%origin_z = init_model%origin_z + oro_min
    1008952
    1009     ENDIF   
    1010    
    1011 !
    1012 !-- In the following, buildings and orography are further preprocessed
    1013 !-- before they are mapped on the LES grid.
    1014 !-- Buildings are mapped on top of the orography by maintaining the roof
    1015 !-- shape of the building. This can be achieved by referencing building on
    1016 !-- top of the maximum terrain height within the area occupied by the
    1017 !-- respective building. As buildings and terrain height are defined PE-wise,
    1018 !-- parallelization of this referencing is required (a building can be
    1019 !-- distributed between different PEs). 
    1020 !-- In a first step, determine the number of buildings with different
    1021 !-- building id on each PE. In a next step, all building ids are gathered
    1022 !-- into one array which is present to all PEs. For each building ID,
    1023 !-- the maximum terrain height occupied by the respective building is
    1024 !-- computed and distributed to each PE. 
    1025 !-- Finally, for each building id and its respective reference orography,
    1026 !-- builidings are mapped on top.   
    1027 !--
    1028 !-- First, pre-set topography flags, bit 1 indicates orography, bit 2
    1029 !-- buildings
    1030 !-- classify the respective surfaces.
     953    ENDIF
     954
     955!
     956!-- In the following, buildings and orography are further preprocessed before they are mapped on the
     957!-- LES grid.
     958!-- Buildings are mapped on top of the orography by maintaining the roof shape of the building. This
     959!-- can be achieved by referencing building on top of the maximum terrain height within the area
     960!-- occupied by the respective building. As buildings and terrain height are defined PE-wise,
     961!-- parallelization of this referencing is required (a building can be distributed between different
     962!-- PEs).
     963!-- In a first step, determine the number of buildings with different building id on each PE. In a
     964!-- next step, all building ids are gathered into one array which is present to all PEs. For each
     965!-- building ID, the maximum terrain height occupied by the respective building is computed and
     966!-- distributed to each PE.
     967!-- Finally, for each building id and its respective reference orography, builidings are mapped on
     968!-- top.
     969!--
     970!-- First, pre-set topography flags, bit 1 indicates orography, bit 2 buildings classify the
     971!-- respective surfaces.
    1031972    topo_3d          = IBSET( topo_3d, 0 )
    1032973    topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
    1033974!
    1034 !-- In order to map topography on PALM grid also in case of ocean simulations,
    1035 !-- pre-calculate an offset value.
     975!-- In order to map topography on PALM grid also in case of ocean simulations, pre-calculate an
     976!-- offset value.
    1036977    ocean_offset = MERGE( zw(0), 0.0_wp, ocean_mode )
    1037978!
    1038 !-- Reference buildings on top of orography. This is not necessary
    1039 !-- if topography is read from ASCII file as no distinction between buildings
    1040 !-- and terrain height can be made. Moreover, this is also not necessary if
    1041 !-- urban-surface and land-surface model are used at the same time.
     979!-- Reference buildings on top of orography. This is not necessary if topography is read from ASCII
     980!-- file as no distinction between buildings and terrain height can be made. Moreover, this is also
     981!-- not necessary if urban-surface and land-surface model are used at the same time.
    1042982    IF ( input_pids_static )  THEN
    1043983
    1044        IF ( buildings_f%from_file )  THEN 
     984       IF ( buildings_f%from_file )  THEN
    1045985          num_buildings_l = 0
    1046986          num_buildings   = 0
    1047987!
    1048 !--       Allocate at least one element for building ids and give it an inital
    1049 !--       negative value that will be overwritten later. This, however, is
    1050 !--       necessary in case there all IDs in the model domain are fill values.
     988!--       Allocate at least one element for building ids and give it an inital negative value that
     989!--       will be overwritten later. This, however, is necessary in case there all IDs in the model
     990!--       domain are fill values.
    1051991          ALLOCATE( build_ids_l(1) )
    1052           build_ids_l = -1 
     992          build_ids_l = -1
    1053993          DO  i = nxl, nxr
    1054994             DO  j = nys, nyn
    1055995                IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
    1056996                   IF ( num_buildings_l(myid) > 0 )  THEN
    1057                       IF ( ANY( building_id_f%var(j,i) ==  build_ids_l ) )   &
    1058                       THEN
     997                      IF ( ANY( building_id_f%var(j,i) ==  build_ids_l ) )  THEN
    1059998                         CYCLE
    1060999                      ELSE
     
    10661005                      DEALLOCATE( build_ids_l )
    10671006                      ALLOCATE( build_ids_l(1:num_buildings_l(myid)) )
    1068                       build_ids_l(1:num_buildings_l(myid)-1) =                 &
    1069                                   build_ids_l_tmp(1:num_buildings_l(myid)-1)
     1007                      build_ids_l(1:num_buildings_l(myid)-1) =                                     &
     1008                                                          build_ids_l_tmp(1:num_buildings_l(myid)-1)
    10701009                      build_ids_l(num_buildings_l(myid)) = building_id_f%var(j,i)
    10711010                      DEALLOCATE( build_ids_l_tmp )
    10721011                   ENDIF
    10731012!
    1074 !--                First occuring building id on PE 
    1075                    ELSE 
     1013!--                First occuring building id on PE
     1014                   ELSE
    10761015                      num_buildings_l(myid) = num_buildings_l(myid) + 1
    10771016                      build_ids_l(1) = building_id_f%var(j,i)
     
    10811020          ENDDO
    10821021!
    1083 !--       Determine number of different building ids for the entire domain 
    1084 #if defined( __parallel ) 
    1085           CALL MPI_ALLREDUCE( num_buildings_l, num_buildings, numprocs,              &
    1086                               MPI_INTEGER, MPI_SUM, comm2d, ierr )
     1022!--       Determine number of different building ids for the entire domain
     1023#if defined( __parallel )
     1024          CALL MPI_ALLREDUCE( num_buildings_l, num_buildings, numprocs, MPI_INTEGER, MPI_SUM,      &
     1025                              comm2d, ierr )
    10871026#else
    10881027          num_buildings = num_buildings_l
    10891028#endif
    10901029!
    1091 !--       Gather all buildings ids on each PEs. 
    1092 !--       First, allocate array encompassing all building ids in model domain. 
     1030!--       Gather all buildings ids on each PEs.
     1031!--       First, allocate array encompassing all building ids in model domain.
    10931032          ALLOCATE( build_ids(1:SUM(num_buildings)) )
    1094 #if defined( __parallel )
    1095 !
    1096 !--       Allocate array for displacements.
    1097 !--       As each PE may has a different number of buildings, so that
    1098 !--       the block sizes send by each PE may not be equal. Hence,
    1099 !--       information about the respective displacement is required, indicating
    1100 !--       the respective adress where each MPI-task writes into the receive
    1101 !--       buffer array 
     1033#if defined( __parallel )
     1034!
     1035!--       Allocate array for displacements.
     1036!--       As each PE may has a different number of buildings, so that the block sizes send by each
     1037!--       PE may not be equal. Hence,  information about the respective displacement is required,
     1038!--       indicating the respective adress where each MPI-task writes into the receive buffer array.
    11021039          ALLOCATE( displace_dum(0:numprocs-1) )
    11031040          displace_dum(0) = 0
     
    11061043          ENDDO
    11071044
    1108           CALL MPI_ALLGATHERV( build_ids_l(1:num_buildings_l(myid)),                 &
    1109                                num_buildings(myid),                                  &
    1110                                MPI_INTEGER,                                          &
    1111                                build_ids,                                            &
    1112                                num_buildings,                                        &
    1113                                displace_dum,                                         &
    1114                                MPI_INTEGER,                                          &
    1115                                comm2d, ierr )   
     1045          CALL MPI_ALLGATHERV( build_ids_l(1:num_buildings_l(myid)), num_buildings(myid),          &
     1046                               MPI_INTEGER, build_ids, num_buildings, displace_dum, MPI_INTEGER,   &
     1047                               comm2d, ierr )
    11161048
    11171049          DEALLOCATE( displace_dum )
     
    11221054
    11231055!
    1124 !--       Note, in parallel mode building ids can occure mutliple times, as
    1125 !--       each PE has send its own ids. Therefore, sort out building ids which
    1126 !--       appear more than one time.
     1056!--       Note, in parallel mode building ids can occure mutliple times, as each PE has send its own
     1057!--       ids. Therefore, sort out building ids which appear more than one time.
    11271058          num_build = 0
    11281059          DO  nr = 1, SIZE(build_ids)
     
    11421073                   build_ids_final(num_build) = build_ids(nr)
    11431074                   DEALLOCATE( build_ids_final_tmp )
    1144                 ENDIF             
     1075                ENDIF
    11451076             ELSE
    11461077                num_build = num_build + 1
     
    11511082
    11521083!
    1153 !--       Determine maximumum terrain height occupied by the respective
    1154 !--       building and temporalily store on oro_max
     1084!--       Determine maximumum terrain height occupied by the respective building and temporalily
     1085!--       store on oro_max.
    11551086          ALLOCATE( oro_max_l(1:SIZE(build_ids_final)) )
    11561087          ALLOCATE( oro_max(1:SIZE(build_ids_final))   )
     
    11581089
    11591090          DO  nr = 1, SIZE(build_ids_final)
    1160              oro_max_l(nr) = MAXVAL(                                           &
    1161                               MERGE( terrain_height_f%var(nys:nyn,nxl:nxr),    &
    1162                                      0.0_wp,                                   &
    1163                                      building_id_f%var(nys:nyn,nxl:nxr) ==     &
    1164                                      build_ids_final(nr) ) )
    1165           ENDDO
    1166    
    1167 #if defined( __parallel )   
    1168           IF ( SIZE(build_ids_final) >= 1 ) THEN
    1169              CALL MPI_ALLREDUCE( oro_max_l, oro_max, SIZE( oro_max ), MPI_REAL,&
    1170                                  MPI_MAX, comm2d, ierr )
     1091             oro_max_l(nr) = MAXVAL( MERGE( terrain_height_f%var(nys:nyn,nxl:nxr),                 &
     1092                                              0.0_wp,                                              &
     1093                                              building_id_f%var(nys:nyn,nxl:nxr) ==                &
     1094                                              build_ids_final(nr) ) )
     1095          ENDDO
     1096
     1097#if defined( __parallel )
     1098          IF ( SIZE(build_ids_final) >= 1 )  THEN
     1099             CALL MPI_ALLREDUCE( oro_max_l, oro_max, SIZE( oro_max ), MPI_REAL, MPI_MAX, comm2d,   &
     1100                                 ierr )
    11711101          ENDIF
    11721102#else
     
    11741104#endif
    11751105!
    1176 !--       Finally, determine discrete grid height of maximum orography occupied
    1177 !--       by a building. Use all-or-nothing approach, i.e. if terrain
    1178 !--       exceeds the scalar level the grid box is fully terrain and the
    1179 !--       maximum terrain is set to the zw level.
    1180 !--       terrain or
     1106!--       Finally, determine discrete grid height of maximum orography occupied by a building. Use
     1107!--       all-or-nothing approach, i.e. if terrain exceeds the scalar level the grid box is fully
     1108!--       terrain and the maximum terrain is set to the zw level.
     1109!--       terrain or
    11811110          oro_max_l = 0.0
    11821111          DO  nr = 1, SIZE(build_ids_final)
    11831112             DO  k = nzb, nzt
    1184                 IF ( zu(k) - ocean_offset <= oro_max(nr) )                     &
    1185                    oro_max_l(nr) = zw(k) - ocean_offset
     1113                IF ( zu(k) - ocean_offset <= oro_max(nr) )  oro_max_l(nr) = zw(k) - ocean_offset
    11861114             ENDDO
    11871115             oro_max(nr) = oro_max_l(nr)
     
    11951123       END IF
    11961124!
    1197 !--    Map orography as well as buildings onto grid. 
     1125!--    Map orography as well as buildings onto grid.
    11981126       DO  i = nxl, nxr
    11991127          DO  j = nys, nyn
     
    12041132                IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
    12051133!
    1206 !--                Determine index where maximum terrain height occupied by
    1207 !--                the respective building height is stored.
    1208                    nr = MINLOC( ABS( build_ids_final -                         &
    1209                                      building_id_f%var(j,i) ), DIM = 1 )
     1134!--                Determine index where maximum terrain height occupied by the respective building
     1135!--                height is stored.
     1136                   nr = MINLOC( ABS( build_ids_final - building_id_f%var(j,i) ), DIM=1 )
    12101137!
    12111138!--                Save grid-indexed oro_max
     
    12151142             DO  k = nzb, nzt
    12161143!
    1217 !--             In a first step, if grid point is below or equal the given
    1218 !--             terrain height, grid point is flagged to be of type natural.
    1219 !--             Please note, in case there is also a building which is lower
    1220 !--             than the vertical grid spacing, initialization of surface
    1221 !--             attributes will not be correct as given surface information
    1222 !--             will not be in accordance to the classified grid points.
     1144!--             In a first step, if grid point is below or equal the given terrain height, grid
     1145!--             point is flagged to be of type natural.
     1146!--             Please note, in case there is also a building which is lower than the vertical grid
     1147!--             spacing, initialization of surface attributes will not be correct as given surface
     1148!--             information will not be in accordance to the classified grid points.
    12231149!--             Hence, in this case, also a building flag.
    12241150                IF ( zu(k) - ocean_offset <= terrain_height_f%var(j,i) )  THEN
     
    12281154                ENDIF
    12291155!
    1230 !--             Set building grid points. Here, only consider 2D buildings. 
    1231 !--             3D buildings require separate treatment. 
     1156!--             Set building grid points. Here, only consider 2D buildings.
     1157!--             3D buildings require separate treatment.
    12321158                IF ( buildings_f%from_file  .AND.  buildings_f%lod == 1 )  THEN
    12331159!
    1234 !--                Fill-up the terrain to the level of maximum orography
    1235 !--                within the building-covered area.
     1160!--                Fill-up the terrain to the level of maximum orography within the building-covered
     1161!--                area.
    12361162                   IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
    12371163!
    1238 !--                   Note, oro_max is always on zw level                   
     1164!--                   Note, oro_max is always on zw level
    12391165                      IF ( zu(k) - ocean_offset < oro_max(nr) )  THEN
    12401166                         topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
    12411167                         topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
    1242                       ELSEIF ( zu(k) - ocean_offset <=                         &
    1243                                oro_max(nr) + buildings_f%var_2d(j,i) )  THEN
     1168                      ELSEIF ( zu(k) - ocean_offset <= oro_max(nr) + buildings_f%var_2d(j,i) )  THEN
    12441169                         topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
    12451170                         topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
     
    12491174             ENDDO
    12501175!
    1251 !--          Special treatment for non grid-resolved buildings. This case,
    1252 !--          the uppermost terrain grid point is flagged as building as well
    1253 !--          well, even though no building exists at all. However, the
    1254 !--          surface element will be identified as urban-surface and the
    1255 !--          input data provided by the drivers is consistent to the surface
    1256 !--          classification. Else, all non grid-resolved buildings would vanish
    1257 !--          and identified as terrain grid points, which, however, won't be
    1258 !--          consistent with the input data.
     1176!--          Special treatment for non grid-resolved buildings. This case, the uppermost terrain
     1177!--          grid point is flagged as building as well, even though no building exists at all.
     1178!--          However, the surface element will be identified as urban-surface and the input data
     1179!--          provided by the drivers is consistent to the surface classification. Else, all non
     1180!--          grid-resolved buildings would vanish and identified as terrain grid points, which,
     1181!--          however, won't be consistent with the input data.
    12591182             IF ( buildings_f%from_file  .AND.  buildings_f%lod == 1 )  THEN
    12601183                IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
     
    12691192             ENDIF
    12701193!
    1271 !--          Map 3D buildings onto terrain height. 
    1272 !--          In case of any slopes, map building on top of maximum terrain
    1273 !--          height covered by the building. In other words, extend
    1274 !--          building down to the respective local terrain-surface height.
     1194!--          Map 3D buildings onto terrain height.
     1195!--          In case of any slopes, map building on top of maximum terrain height covered by the
     1196!--          building. In other words, extend building down to the respective local terrain-surface
     1197!--          height.
    12751198             IF ( buildings_f%from_file  .AND.  buildings_f%lod == 2 )  THEN
    12761199                IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
    12771200!
    1278 !--                Extend building down to the terrain surface, i.e. fill-up
    1279 !--                surface irregularities below a building. Note, oro_max
    1280 !--                is already a discrete height according to the all-or-nothing
    1281 !--                approach, i.e. grid box is either topography or atmosphere,
     1201!--                Extend building down to the terrain surface, i.e. fill-up surface irregularities
     1202!--                below a building. Note, oro_max is already a discrete height according to the
     1203!--                all-or-nothing approach, i.e. grid box is either topography or atmosphere,
    12821204!--                terrain top is defined at upper bound of the grid box.
    1283 !--                Hence, check for zw in this case.
    1284 !--                Note, do this only for buildings which are surface mounted,
    1285 !--                i.e. building types 1-6. Below bridges, which are represented
    1286 !--                exclusively by building type 7, terrain shape should be
    1287 !--                maintained.
     1205!--                Hence, check for zw in this case.
     1206!--                Note, do this only for buildings which are surface mounted, i.e. building types
     1207!--                1-6. Below bridges, which are represented exclusively by building type 7, terrain
     1208!--                shape should be maintained.
    12881209                   IF ( building_type_f%from_file )  THEN
    12891210                      IF ( building_type_f%var(j,i) /= 7 )  THEN
    1290                          DO k = topo_top_index + 1, nzt + 1     
     1211                         DO k = topo_top_index + 1, nzt + 1
    12911212                            IF ( zu(k) - ocean_offset <= oro_max(nr) )  THEN
    12921213                               topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
    12931214                               topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
    12941215                            ENDIF
    1295                          ENDDO       
    1296 !                     
    1297 !--                      After surface irregularities are smoothen, determine
    1298 !--                      lower start index where building starts.
     1216                         ENDDO
     1217!
     1218!--                      After surface irregularities are smoothen, determine lower start index
     1219!--                      where building starts.
    12991220                         DO  k = nzb, nzt
    1300                             IF ( zu(k) - ocean_offset <= oro_max(nr) )         &
    1301                                topo_top_index = k
     1221                            IF ( zu(k) - ocean_offset <= oro_max(nr) )  topo_top_index = k
    13021222                         ENDDO
    13031223                      ENDIF
     
    13201240       ENDDO
    13211241!
    1322 !--    Horizontal exchange the oro_max array, which is required to for
    1323 !--    initialization of building-surface properties.
     1242!--    Horizontal exchange the oro_max array, which is required to for initialization of
     1243!--    building-surface properties.
    13241244       IF ( ALLOCATED( buildings_f%oro_max ) )  THEN
    13251245          CALL exchange_horiz_2d( buildings_f%oro_max(:,:) )
     
    13311251       IF ( ALLOCATED( build_ids_final ) )  DEALLOCATE( build_ids_final )
    13321252!
    1333 !-- Topography input via ASCII format. 
     1253!-- Topography input via ASCII format.
    13341254    ELSE
    13351255       ocean_offset     = MERGE( zw(0), 0.0_wp, ocean_mode )
    13361256!
    1337 !--    Initialize topography bit 0 (indicates obstacle) everywhere to zero
    1338 !--    and clear all grid points at nzb, where alway a surface is defined.
    1339 !--    Further, set also bit 1 (indicates terrain) at nzb, which is further
    1340 !--    used for masked data output and further processing. Note, in the
    1341 !--    ASCII case no distinction is made between buildings and terrain,
    1342 !--    so that setting of bit 1 and 2 at the same time has no effect.
     1257!--    Initialize topography bit 0 (indicates obstacle) everywhere to zero and clear all grid points
     1258!--    at nzb, where alway a surface is defined.
     1259!--    Further, set also bit 1 (indicates terrain) at nzb, which is further used for masked data
     1260!--    output and further processing. Note, in the ASCII case no distinction is made between
     1261!--    buildings and terrain,  so that setting of bit 1 and 2 at the same time has no effect.
    13431262       topo_3d          = IBSET( topo_3d, 0 )
    13441263       topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
     
    13481267             DO  k = nzb, nzt
    13491268!
    1350 !--             Flag topography for all grid points which are below
    1351 !--             the local topography height.
    1352 !--             Note, each topography is flagged as building (bit 2) as well as
    1353 !--             terrain (bit 1) in order to employ urban-surface as well as
    1354 !--             land-surface model.
     1269!--             Flag topography for all grid points which are below the local topography height.
     1270!--             Note, each topography is flagged as building (bit 2) as well as terrain (bit 1) in
     1271!--             order to employ urban-surface as well as land-surface model.
    13551272                IF ( zu(k) - ocean_offset <= buildings_f%var_2d(j,i) )  THEN
    13561273                   topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
     
    13721289    IF ( .NOT. bc_lr_cyc )  THEN
    13731290       IF ( nxl == 0  )  topo_3d(:,:,-1)   = topo_3d(:,:,0)
    1374        IF ( nxr == nx )  topo_3d(:,:,nx+1) = topo_3d(:,:,nx)         
     1291       IF ( nxr == nx )  topo_3d(:,:,nx+1) = topo_3d(:,:,nx)
    13751292    ENDIF
    13761293
     
    13791296
    13801297! Description:
    1381 ! -----------------------------------------------------------------------------!
    1382 !> Filter topography, i.e. fill holes resolved by only one grid point. 
    1383 !> Such holes are suspected to lead to velocity blow-ups as continuity
    1384 !> equation on discrete grid cannot be fulfilled in such case.
    1385 !------------------------------------------------------------------------------!
     1298! -------------------------------------------------------------------------------------------------!
     1299!> Filter topography, i.e. fill holes resolved by only one grid point.
     1300!> Such holes are suspected to lead to velocity blow-ups as continuity equation on discrete grid
     1301!> cannot be fulfilled in such case.
     1302!--------------------------------------------------------------------------------------------------!
    13861303 SUBROUTINE filter_topography( topo_3d )
    13871304
    1388     USE control_parameters,                                                    &
     1305    USE control_parameters,                                                                        &
    13891306        ONLY:  bc_lr_cyc, bc_ns_cyc, message_string
    13901307
    1391     USE exchange_horiz_mod,                                                    &
     1308    USE exchange_horiz_mod,                                                                        &
    13921309        ONLY:  exchange_horiz_int, exchange_horiz_2d_byte, exchange_horiz_2d_int
    13931310
    1394     USE indices,                                                               &
     1311    USE indices,                                                                                   &
    13951312        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nzt
    13961313
    1397     USE netcdf_data_input_mod,                                                 &
    1398         ONLY:  building_id_f, building_type_f 
     1314    USE netcdf_data_input_mod,                                                                     &
     1315        ONLY:  building_id_f, building_type_f
    13991316
    14001317    USE  pegrid
    14011318
    14021319    IMPLICIT NONE
    1403 
    1404     LOGICAL      ::  filled = .FALSE. !< flag indicating if holes were filled
    14051320
    14061321    INTEGER(iwp) ::  i          !< running index along x-direction
    14071322    INTEGER(iwp) ::  j          !< running index along y-direction
    14081323    INTEGER(iwp) ::  k          !< running index along z-direction
    1409     INTEGER(iwp) ::  num_hole   !< number of holes (in topography) resolved by only one grid point 
    1410     INTEGER(iwp) ::  num_hole_l !< number of holes (in topography) resolved by only one grid point on local PE     
     1324    INTEGER(iwp) ::  num_hole   !< number of holes (in topography) resolved by only one grid point
     1325    INTEGER(iwp) ::  num_hole_l !< number of holes (in topography) resolved by only one grid point on local PE
    14111326    INTEGER(iwp) ::  num_wall   !< number of surrounding vertical walls for a single grid point
    14121327
    14131328    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE            ::  topo_tmp          !< temporary 3D-topography used to fill holes
    1414     INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo_3d           !< 3D-topography array merging buildings and orography
    1415 !
    1416 !-- Before checking for holes, set lateral boundary conditions for
     1329    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo_3d           !< 3D-topography array merging buildings and
     1330                                                                                 !< orography
     1331
     1332    LOGICAL      ::  filled = .FALSE. !< flag indicating if holes were filled
     1333
     1334!
     1335!-- Before checking for holes, set lateral boundary conditions for
    14171336!-- topography. After hole-filling, boundary conditions must be set again.
    1418 !-- Several iterations are performed, in order to fill holes which might 
     1337!-- Several iterations are performed, in order to fill holes which might
    14191338!-- emerge by the filling-algorithm itself.
    14201339    ALLOCATE( topo_tmp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     
    14221341
    14231342    num_hole = 99999
    1424     DO WHILE ( num_hole > 0 )       
    1425 
    1426        num_hole = 0   
     1343    DO WHILE ( num_hole > 0 )
     1344
     1345       num_hole = 0
    14271346       CALL exchange_horiz_int( topo_3d, nys, nyn, nxl, nxr, nzt, nbgp )
    14281347!
    1429 !--    Exchange also building ID and type. Note, building_type is an one-byte
    1430 !--    variable.
    1431        IF ( building_id_f%from_file )                                          &
     1348!--    Exchange also building ID and type. Note, building_type is an one-byte variable.
     1349       IF ( building_id_f%from_file )                                                              &
    14321350          CALL exchange_horiz_2d_int( building_id_f%var, nys, nyn, nxl, nxr, nbgp )
    1433        IF ( building_type_f%from_file )                                        &
     1351       IF ( building_type_f%from_file )                                                            &
    14341352          CALL exchange_horiz_2d_byte( building_type_f%var, nys, nyn, nxl, nxr, nbgp )
    14351353
    14361354       topo_tmp = topo_3d
    14371355!
    1438 !--    In case of non-cyclic lateral boundaries, assume lateral boundary to be
    1439 !--    a solid wall. Thus, intermediate spaces of one grid point between
    1440 !--    boundary and some topographic structure will be filled.           
     1356!--    In case of non-cyclic lateral boundaries, assume lateral boundary to be a solid wall. Thus,
     1357!--    intermediate spaces of one grid point between boundary and some topographic structure will be
     1358!--    filled.
    14411359       IF ( .NOT. bc_ns_cyc )  THEN
    14421360          IF ( nys == 0  )  topo_tmp(:,-1,:)   = IBCLR( topo_tmp(:,0,:),  0 )
     
    14461364       IF ( .NOT. bc_lr_cyc )  THEN
    14471365          IF ( nxl == 0  )  topo_tmp(:,:,-1)   = IBCLR( topo_tmp(:,:,0),  0 )
    1448           IF ( nxr == nx )  topo_tmp(:,:,nx+1) = IBCLR( topo_tmp(:,:,nx), 0 )         
     1366          IF ( nxr == nx )  topo_tmp(:,:,nx+1) = IBCLR( topo_tmp(:,:,nx), 0 )
    14491367       ENDIF
    14501368
     
    14551373                IF ( BTEST( topo_tmp(k,j,i), 0 ) )  THEN
    14561374                   num_wall = 0
    1457                    IF ( .NOT. BTEST( topo_tmp(k,j-1,i), 0 ) )                  &
    1458                       num_wall = num_wall + 1
    1459                    IF ( .NOT. BTEST( topo_tmp(k,j+1,i), 0 ) )                  &
    1460                       num_wall = num_wall + 1
    1461                    IF ( .NOT. BTEST( topo_tmp(k,j,i-1), 0 ) )                  &
    1462                       num_wall = num_wall + 1
    1463                    IF ( .NOT. BTEST( topo_tmp(k,j,i+1), 0 ) )                  &
    1464                       num_wall = num_wall + 1
    1465                    IF ( .NOT. BTEST( topo_tmp(k-1,j,i), 0 ) )                  &
    1466                       num_wall = num_wall + 1   
    1467                    IF ( .NOT. BTEST( topo_tmp(k+1,j,i), 0 ) )                  &
    1468                       num_wall = num_wall + 1
     1375                   IF ( .NOT. BTEST( topo_tmp(k,j-1,i), 0 ) )  num_wall = num_wall + 1
     1376                   IF ( .NOT. BTEST( topo_tmp(k,j+1,i), 0 ) )  num_wall = num_wall + 1
     1377                   IF ( .NOT. BTEST( topo_tmp(k,j,i-1), 0 ) )  num_wall = num_wall + 1
     1378                   IF ( .NOT. BTEST( topo_tmp(k,j,i+1), 0 ) )  num_wall = num_wall + 1
     1379                   IF ( .NOT. BTEST( topo_tmp(k-1,j,i), 0 ) )  num_wall = num_wall + 1
     1380                   IF ( .NOT. BTEST( topo_tmp(k+1,j,i), 0 ) )  num_wall = num_wall + 1
    14691381
    14701382                   IF ( num_wall >= 4 )  THEN
    14711383                      num_hole_l     = num_hole_l + 1
    14721384!
    1473 !--                   Clear flag 0 and set special flag ( bit 4) to indicate
    1474 !--                   that new topography point is a result of filtering process.
     1385!--                   Clear flag 0 and set special flag ( bit 4) to indicate that new topography
     1386!--                   point is a result of filtering process.
    14751387                      topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
    14761388                      topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 4 )
    14771389!
    1478 !--                   If filled grid point is occupied by a building, classify
    1479 !--                   it as building grid point.
     1390!--                   If filled grid point is occupied by a building, classify it as building grid
     1391!--                   point.
    14801392                      IF ( building_type_f%from_file )  THEN
    1481                          IF ( building_type_f%var(j,i)   /=                    & 
    1482                               building_type_f%fill            .OR.             &       
    1483                               building_type_f%var(j+1,i) /=                    & 
    1484                               building_type_f%fill            .OR.             &               
    1485                               building_type_f%var(j-1,i) /=                    &               
    1486                               building_type_f%fill            .OR.             &               
    1487                               building_type_f%var(j,i+1) /=                    &               
    1488                               building_type_f%fill            .OR.             &               
    1489                               building_type_f%var(j,i-1) /=                    &               
    1490                               building_type_f%fill )  THEN
     1393                         IF ( building_type_f%var(j,i)   /=  building_type_f%fill  .OR.            &
     1394                              building_type_f%var(j+1,i) /=  building_type_f%fill  .OR.            &
     1395                              building_type_f%var(j-1,i) /=  building_type_f%fill  .OR.            &
     1396                              building_type_f%var(j,i+1) /=  building_type_f%fill  .OR.            &
     1397                              building_type_f%var(j,i-1) /=  building_type_f%fill )  THEN
    14911398!
    14921399!--                         Set flag indicating building surfaces
    14931400                            topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
    14941401!
    1495 !--                         Set building_type and ID at this position if not
    1496 !--                         already set. This is required for proper
    1497 !--                         initialization of urban-surface energy balance
     1402!--                         Set building_type and ID at this position if not already set. This is
     1403!--                         required for proper initialization of urban-surface energy balance
    14981404!--                         solver.
    1499                             IF ( building_type_f%var(j,i) ==                   &
    1500                                  building_type_f%fill )  THEN
    1501 
    1502                                IF ( building_type_f%var(j+1,i) /=              &
    1503                                     building_type_f%fill )  THEN
    1504                                   building_type_f%var(j,i) =                   &
    1505                                                     building_type_f%var(j+1,i)
    1506                                   building_id_f%var(j,i) =                     &
    1507                                                     building_id_f%var(j+1,i)
    1508