Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

files re-formatted to follow the PALM coding standard

File:
1 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 )
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |