Changeset 4803


Ignore:
Timestamp:
Nov 30, 2020 4:57:54 PM (4 years ago)
Author:
raasch
Message:

file re-formatted to follow the PALM coding standard

File:
1 edited

Legend:

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

    r4770 r4803  
    11!> @file plant_canopy_model_mod.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
    1817! Copyright 2017-2020 Institute of Computer Science of the
    1918!                     Czech Academy of Sciences, Prague
    20 !------------------------------------------------------------------------------!
     19!--------------------------------------------------------------------------------------------------!
    2120!
    2221! Current revisions:
     
    2726! -----------------
    2827! $Id$
    29 ! Consider basal-area density as an additional sink for momentum in the
    30 ! prognostic equations for momentum and SGS TKE
     28! file re-formatted to follow the PALM coding standard
     29!
     30! 4770 2020-11-03 14:04:53Z suehring
     31! Consider basal-area density as an additional sink for momentum in the prognostic equations for
     32! momentum and SGS TKE
    3133!
    3234! 4768 2020-11-02 19:11:23Z suehring
     
    3537! 4671 2020-09-09 20:27:58Z pavelkrc
    3638! Implementation of downward facing USM and LSM surfaces
    37 ! 
     39!
    3840! 4535 2020-05-15 12:07:23Z raasch
    3941! bugfix for restart data format query
    40 ! 
     42!
    4143! 4525 2020-05-10 17:05:07Z raasch
    4244! bugfix for reading/writing pcm_...rate_av with MPI-IO
    43 ! 
     45!
    4446! 4517 2020-05-03 14:29:30Z raasch
    4547! added restart with MPI-IO for reading local arrays
    46 ! 
     48!
    4749! 4515 2020-04-30 16:37:18Z suehring
    4850! Rename error number again since this was recently given in -r 4511
    49 ! 
     51!
    5052! 4514 2020-04-30 16:29:59Z suehring
    51 ! - Bugfix in output of pcm_heatrate_av in a restart run. In order to fix this,
    52 !   pch_index is now output for a restart run. Therefore, define global restart
    53 !   routines.
    54 ! - Error message number renamed and check for PA0505 revised in order to also
    55 !   consider natural surfaces with plant-canopy.
    56 !
     53! - Bugfix in output of pcm_heatrate_av in a restart run. In order to fix this, pch_index is now
     54!   output for a restart run. Therefore, define global restart routines.
     55! - Error message number renamed and check for PA0505 revised in order to also consider natural
     56!   surfaces with plant-canopy.
     57!
    5758! 4495 2020-04-13 20:11:20Z raasch
    5859! restart data handling with MPI-IO added
     
    6364! (salim) removed the error message PA0672 to consider PC 3d data via ascii file
    6465!
    65 ! 4392 2020-01-31 16:14:57Z pavelkrc (resler) 
    66 ! Make pcm_heatrate_av, pcm_latentrate_av public to allow calculation
    67 ! of averaged Bowen ratio in the user procedure
     66! 4392 2020-01-31 16:14:57Z pavelkrc (resler)
     67! Make pcm_heatrate_av, pcm_latentrate_av public to allow calculation of averaged Bowen ratio in the
     68! user procedure
    6869!
    6970! 4381 2020-01-20 13:51:46Z suehring
    7071! Give error message 313 only once
    71 ! 
     72!
    7273! 4363 2020-01-07 18:11:28Z suehring
    7374! Fix for last commit
    74 ! 
     75!
    7576! 4362 2020-01-07 17:15:02Z suehring
    76 ! Input of plant canopy variables from static driver moved to plant-canopy
    77 ! model
    78 !
     77! Input of plant canopy variables from static driver moved to plant-canopy model
     78!
    7979! 4361 2020-01-07 12:22:38Z suehring
    8080! - Remove unused arrays in pmc_rrd_local
    8181! - Remove one exchange of ghost points
    82 ! 
     82!
    8383! 4360 2020-01-07 11:25:50Z suehring
    8484! - Bugfix, read restart data for time-averaged pcm output quantities
    8585! - Output of plant-canopy quantities will fill values
    86 ! 
     86!
    8787! 4356 2019-12-20 17:09:33Z suehring
    88 ! Correct single message call, local check must be given by the respective
    89 ! mpi rank.
    90 !
     88! Correct single message call, local check must be given by the respective mpi rank.
     89!
    9190! 4346 2019-12-18 11:55:56Z motisi
    92 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    93 ! topography information used in wall_flags_static_0
    94 ! 
     91! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     92! information used in wall_flags_static_0
     93!
    9594! 4342 2019-12-16 13:49:14Z Giersch
    96 ! Use statements moved to module level, ocean dependency removed, redundant
    97 ! variables removed
    98 !
     95! Use statements moved to module level, ocean dependency removed, redundant variables removed
     96!
    9997! 4341 2019-12-16 10:43:49Z motisi
    100 ! - Unification of variable names: pc_-variables now pcm_-variables
    101 !   (pc_latent_rate, pc_heating_rate, pc_transpiration_rate)
     98! - Unification of variable names: pc_-variables now pcm_-variables (pc_latent_rate,
     99!                                 pc_heating_rate, pc_transpiration_rate)
    102100! - Removal of pcm_bowenratio output
    103101! - Renamed canopy-mode 'block' to 'homogeneous'
     
    106104! - Replacement of k_wall by topo_top_ind
    107105! - Removal of Else-Statement in tendency-calculation
    108 ! 
     106!
    109107! 4335 2019-12-12 16:39:05Z suehring
    110108! Fix for LAD at building edges also implemented in vector branch.
    111 ! 
     109!
    112110! 4331 2019-12-10 18:25:02Z suehring
    113111! Typo corrected
    114 ! 
     112!
    115113! 4329 2019-12-10 15:46:36Z motisi
    116114! Renamed wall_flags_0 to wall_flags_static_0
    117 ! 
     115!
    118116! 4314 2019-11-29 10:29:20Z suehring
    119 ! - Bugfix, plant canopy was still considered at building edges on for the u-
    120 !   and v-component.
    121 ! - Relax restriction of LAD on building tops. LAD is only omitted at
    122 !   locations where building grid points emerged artificially by the
    123 !   topography filtering.
    124 !
     117! - Bugfix, plant canopy was still considered at building edges on for the u- and v-component.
     118! - Relax restriction of LAD on building tops. LAD is only omitted at locations where building grid
     119!   points emerged artificially by the topography filtering.
     120!
    125121! 4309 2019-11-26 18:49:59Z suehring
    126122! Typo
    127 ! 
     123!
    128124! 4302 2019-11-22 13:15:56Z suehring
    129125! Omit tall canopy mapped on top of buildings
    130 ! 
     126!
    131127! 4279 2019-10-29 08:48:17Z scharf
    132128! unused variables removed
    133 ! 
     129!
    134130! 4258 2019-10-07 13:29:08Z scharf
    135131! changed check for static driver and fixed bugs in initialization and header
    136 ! 
     132!
    137133! 4258 2019-10-07 13:29:08Z suehring
    138134! Check if any LAD is prescribed when plant-canopy model is applied.
    139 ! 
     135!
    140136! 4226 2019-09-10 17:03:24Z suehring
    141137! Bugfix, missing initialization of heating rate
    142 ! 
     138!
    143139! 4221 2019-09-09 08:50:35Z suehring
    144140! Further bugfix in 3d data output for plant canopy
    145 ! 
     141!
    146142! 4216 2019-09-04 09:09:03Z suehring
    147143! Bugfixes in 3d data output
    148 ! 
     144!
    149145! 4205 2019-08-30 13:25:00Z suehring
    150146! Missing working precision + bugfix in calculation of wind speed
    151 ! 
     147!
    152148! 4188 2019-08-26 14:15:47Z suehring
    153149! Minor adjustment in error number
    154 ! 
     150!
    155151! 4187 2019-08-26 12:43:15Z suehring
    156152! Give specific error numbers instead of PA0999
    157 ! 
     153!
    158154! 4182 2019-08-22 15:20:23Z scharf
    159155! Corrected "Former revisions" section
    160 ! 
     156!
    161157! 4168 2019-08-16 13:50:17Z suehring
    162158! Replace function get_topography_top_index by topo_top_ind
    163 ! 
     159!
    164160! 4127 2019-07-30 14:47:10Z suehring
    165 ! Output of 3D plant canopy variables changed. It is now relative to the local
    166 ! terrain rather than located at the acutal vertical level in the model. This
    167 ! way, the vertical dimension of the output can be significantly reduced.
    168 ! (merge from branch resler) 
    169 ! 
     161! Output of 3D plant canopy variables changed. It is now relative to the local terrain rather than
     162! located at the acutal vertical level in the model. This way, the vertical dimension of the output
     163! can be significantly reduced.
     164! (merge from branch resler)
     165!
    170166! 3885 2019-04-11 11:29:34Z kanani
    171 ! Changes related to global restructuring of location messages and introduction
    172 ! of additional debug messages
    173 ! 
     167! Changes related to global restructuring of location messages and introduction of additional debug
     168! messages
     169!
    174170! 3864 2019-04-05 09:01:56Z monakurppa
    175171! unsed variables removed
    176 ! 
     172!
    177173! 3745 2019-02-15 18:57:56Z suehring
    178 ! Bugfix in transpiration, floating invalid when temperature
    179 ! becomes > 40 degrees
    180 !
     174! Bugfix in transpiration, floating invalid when temperature becomes > 40 degrees
     175!
    181176! 3744 2019-02-15 18:38:58Z suehring
    182177! Some interface calls moved to module_interface + cleanup
    183 ! 
     178!
    184179! 3655 2019-01-07 16:51:22Z knoop
    185180! unused variables removed
     
    190185! Description:
    191186! ------------
    192 !> 1) Initialization of the canopy model, e.g. construction of leaf area density
    193 !> profile (subroutine pcm_init).
    194 !> 2) Calculation of sinks and sources of momentum, heat and scalar concentration
    195 !> due to canopy elements (subroutine pcm_tendency).
     187!> 1) Initialization of the canopy model, e.g. construction of leaf area density profile
     188!>   (subroutine pcm_init).
     189!> 2) Calculation of sinks and sources of momentum, heat and scalar concentration due to canopy
     190!>   elements (subroutine pcm_tendency).
    196191!
    197192! @todo - precalculate constant terms in pcm_calc_transpiration_rate
    198193! @todo - unify variable names (pcm_, pc_, ...)
    199194! @todo - get rid-off dependency on radiation model
    200 !------------------------------------------------------------------------------!
     195!--------------------------------------------------------------------------------------------------!
    201196 MODULE plant_canopy_model_mod
    202197
    203     USE arrays_3d,                                                             &
     198    USE arrays_3d,                                                                                 &
    204199        ONLY:  dzu, dzw, e, exner, hyp, pt, q, s, tend, u, v, w, zu, zw
    205200
    206     USE basic_constants_and_equations_mod,                                     &
     201    USE basic_constants_and_equations_mod,                                                         &
    207202        ONLY:  c_p, degc_to_k, l_v, lv_d_cp, r_d, rd_d_rv
    208203
    209     USE bulk_cloud_model_mod,                                                  &
     204    USE bulk_cloud_model_mod,                                                                      &
    210205        ONLY: bulk_cloud_model, microphysics_seifert
    211206
    212     USE control_parameters,                                                    &
    213         ONLY: average_count_3d,                                                &
    214               coupling_char,                                                   &
    215               debug_output,                                                    &
    216               dt_3d,                                                           &
    217               dz,                                                              &
    218               humidity,                                                        &
    219               land_surface,                                                    &
    220               length,                                                          &
    221               message_string,                                                  &
    222               ocean_mode,                                                      &
    223               passive_scalar,                                                  &
    224               plant_canopy,                                                    &
    225               restart_data_format_output,                                      &
    226               restart_string,                                                  &
     207    USE control_parameters,                                                                        &
     208        ONLY: average_count_3d,                                                                    &
     209              coupling_char,                                                                       &
     210              debug_output,                                                                        &
     211              dt_3d,                                                                               &
     212              dz,                                                                                  &
     213              humidity,                                                                            &
     214              land_surface,                                                                        &
     215              length,                                                                              &
     216              message_string,                                                                      &
     217              ocean_mode,                                                                          &
     218              passive_scalar,                                                                      &
     219              plant_canopy,                                                                        &
     220              restart_data_format_output,                                                          &
     221              restart_string,                                                                      &
    227222              urban_surface
    228223
    229     USE grid_variables,                                                        &
     224    USE grid_variables,                                                                            &
    230225        ONLY:  dx, dy
    231226
    232     USE indices,                                                               &
    233         ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
    234                nz, nzb, nzt, topo_top_ind, wall_flags_total_0
     227    USE indices,                                                                                   &
     228        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, nz, nzb, nzt,         &
     229               topo_top_ind, wall_flags_total_0
    235230
    236231    USE kinds
    237232
    238     USE netcdf_data_input_mod,                                                 &
    239         ONLY:  input_pids_static,                                              &
    240                char_fill,                                                      &
    241                check_existence,                                                &
    242                close_input_file,                                               &
    243                get_attribute,                                                  &
    244                get_dimension_length,                                           &
    245                get_variable,                                                   &
    246                inquire_num_variables,                                          &
    247                inquire_variable_names,                                         &
    248                input_file_static,                                              &
    249                num_var_pids,                                                   &
    250                open_read_file,                                                 &
    251                pids_id,                                                        &
    252                real_3d,                                                        &
     233    USE netcdf_data_input_mod,                                                                     &
     234        ONLY:  char_fill,                                                                          &
     235               check_existence,                                                                    &
     236               close_input_file,                                                                   &
     237               get_attribute,                                                                      &
     238               get_dimension_length,                                                               &
     239               get_variable,                                                                       &
     240               input_file_static,                                                                  &
     241               input_pids_static,                                                                  &
     242               inquire_num_variables,                                                              &
     243               inquire_variable_names,                                                             &
     244               num_var_pids,                                                                       &
     245               open_read_file,                                                                     &
     246               pids_id,                                                                            &
     247               real_3d,                                                                            &
    253248               vars_pids
    254249
     
    260255               wrd_mpi_io
    261256
    262     USE surface_mod,                                                           &
     257    USE surface_mod,                                                                               &
    263258        ONLY: surf_def_h, surf_lsm_h, surf_usm_h
    264259
     
    266261    IMPLICIT NONE
    267262
    268     CHARACTER (LEN=30)   ::  canopy_mode = 'homogeneous'           !< canopy coverage
    269     LOGICAL              ::  plant_canopy_transpiration = .FALSE.  !< flag to switch calculation of transpiration and corresponding latent heat
    270                                                                    !< for resolved plant canopy inside radiation model
    271                                                                    !< (calls subroutine pcm_calc_transpiration_rate from module plant_canopy_mod)
    272 
     263    CHARACTER (LEN=30) ::  canopy_mode = 'homogeneous'           !< canopy coverage
    273264    INTEGER(iwp) ::  pch_index = 0                                 !< plant canopy height/top index
     265
    274266    INTEGER(iwp) ::  lad_vertical_gradient_level_ind(10) = -9999   !< lad-profile levels (index)
    275267
    276268    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pch_index_ji     !< local plant canopy top
    277269
    278     LOGICAL ::  calc_beta_lad_profile = .FALSE.   !< switch for calc. of lad from beta func.
    279 
    280     REAL(wp) ::  alpha_lad = 9999999.9_wp         !< coefficient for lad calculation
    281     REAL(wp) ::  beta_lad = 9999999.9_wp          !< coefficient for lad calculation
    282     REAL(wp) ::  canopy_drag_coeff = 0.0_wp       !< canopy drag coefficient (parameter)
    283     REAL(wp) ::  cthf = 0.0_wp                    !< canopy top heat flux
    284     REAL(wp) ::  dt_plant_canopy = 0.0_wp         !< timestep account. for canopy drag
    285     REAL(wp) ::  ext_coef = 0.6_wp                !< extinction coefficient
    286     REAL(wp) ::  lad_surface = 0.0_wp             !< lad surface value
    287     REAL(wp) ::  lai_beta = 0.0_wp                !< leaf area index (lai) for lad calc.
    288     REAL(wp) ::  leaf_scalar_exch_coeff = 0.0_wp  !< canopy scalar exchange coeff.
    289     REAL(wp) ::  leaf_surface_conc = 0.0_wp       !< leaf surface concentration
    290 
     270    LOGICAL ::  calc_beta_lad_profile = .FALSE.       !< switch for calc. of lad from beta func.
     271    LOGICAL ::  plant_canopy_transpiration = .FALSE.  !< flag to switch calculation of transpiration and corresponding latent heat
     272                                                      !< for resolved plant canopy inside radiation model
     273                                                      !< (calls subroutine pcm_calc_transpiration_rate from module plant_canopy_mod)
     274
     275    REAL(wp) ::  alpha_lad = 9999999.9_wp                        !< coefficient for lad calculation
     276    REAL(wp) ::  beta_lad = 9999999.9_wp                         !< coefficient for lad calculation
     277    REAL(wp) ::  canopy_drag_coeff = 0.0_wp                      !< canopy drag coefficient (parameter)
     278    REAL(wp) ::  cthf = 0.0_wp                                   !< canopy top heat flux
     279    REAL(wp) ::  dt_plant_canopy = 0.0_wp                        !< timestep account. for canopy drag
     280    REAL(wp) ::  ext_coef = 0.6_wp                               !< extinction coefficient
     281    REAL(wp) ::  lad_surface = 0.0_wp                            !< lad surface value
     282    REAL(wp) ::  lad_type_coef(0:10) = 1.0_wp                    !< multiplicative coeficients for particular types
     283                                                                 !< of plant canopy (e.g. deciduous tree during winter)
    291284    REAL(wp) ::  lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
    292285    REAL(wp) ::  lad_vertical_gradient_level(10) = -9999999.9_wp !< lad-prof. levels (in m)
    293 
    294     REAL(wp) ::  lad_type_coef(0:10) = 1.0_wp   !< multiplicative coeficients for particular types
    295                                                 !< of plant canopy (e.g. deciduous tree during winter)
     286    REAL(wp) ::  lai_beta = 0.0_wp                               !< leaf area index (lai) for lad calc.
     287    REAL(wp) ::  leaf_scalar_exch_coeff = 0.0_wp                 !< canopy scalar exchange coeff.
     288    REAL(wp) ::  leaf_surface_conc = 0.0_wp                      !< leaf surface concentration
    296289
    297290    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !< leaf area density
     
    302295    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s                    !< lad on scalar-grid
    303296    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_heating_rate         !< plant canopy heating rate
     297    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_heatrate_av          !< array for averaging plant canopy sensible heating rate
     298    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_latent_rate          !< plant canopy latent heating rate
     299    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_latentrate_av        !< array for averaging plant canopy latent heating rate
    304300    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_transpiration_rate   !< plant canopy transpiration rate
    305     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_latent_rate          !< plant canopy latent heating rate
    306 
    307     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_heatrate_av          !< array for averaging plant canopy sensible heating rate
    308     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_latentrate_av        !< array for averaging plant canopy latent heating rate
    309301    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_transpirationrate_av !< array for averaging plant canopy transpiration rate
    310302
     
    316308
    317309    PRIVATE
    318    
     310
    319311!
    320312!-- Public functions
    321     PUBLIC pcm_calc_transpiration_rate,                                       &
    322            pcm_check_data_output,                                             &
    323            pcm_check_parameters,                                              &
    324            pcm_3d_data_averaging,                                             &
    325            pcm_data_output_3d,                                                &
    326            pcm_define_netcdf_grid,                                            &
    327            pcm_header,                                                        &
    328            pcm_init,                                                          &
    329            pcm_parin,                                                         &
    330            pcm_rrd_global,                                                    &
    331            pcm_rrd_local,                                                     &
    332            pcm_tendency,                                                      &
    333            pcm_wrd_global,                                                    &
     313    PUBLIC pcm_calc_transpiration_rate,                                                            &
     314           pcm_check_data_output,                                                                  &
     315           pcm_check_parameters,                                                                   &
     316           pcm_3d_data_averaging,                                                                  &
     317           pcm_data_output_3d,                                                                     &
     318           pcm_define_netcdf_grid,                                                                 &
     319           pcm_header,                                                                             &
     320           pcm_init,                                                                               &
     321           pcm_parin,                                                                              &
     322           pcm_rrd_global,                                                                         &
     323           pcm_rrd_local,                                                                          &
     324           pcm_tendency,                                                                           &
     325           pcm_wrd_global,                                                                         &
    334326           pcm_wrd_local
    335327
    336328!
    337329!-- Public variables and constants
    338     PUBLIC canopy_drag_coeff, pcm_heating_rate, pcm_transpiration_rate,       &
    339            pcm_latent_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s,   &
    340            pch_index, plant_canopy_transpiration,                             &
     330    PUBLIC canopy_drag_coeff, pcm_heating_rate, pcm_transpiration_rate, pcm_latent_rate,           &
     331           canopy_mode, cthf, dt_plant_canopy, lad, lad_s, pch_index, plant_canopy_transpiration,  &
    341332           pcm_heatrate_av, pcm_latentrate_av
    342333
     
    348339       MODULE PROCEDURE pcm_check_data_output
    349340    END INTERFACE pcm_check_data_output
    350    
     341
    351342    INTERFACE pcm_check_parameters
    352343       MODULE PROCEDURE pcm_check_parameters
     
    364355       MODULE PROCEDURE pcm_define_netcdf_grid
    365356    END INTERFACE pcm_define_netcdf_grid
    366    
     357
    367358     INTERFACE pcm_header
    368359       MODULE PROCEDURE pcm_header
    369     END INTERFACE pcm_header       
    370    
     360    END INTERFACE pcm_header
     361
    371362    INTERFACE pcm_init
    372363       MODULE PROCEDURE pcm_init
     
    406397
    407398 CONTAINS
    408  
    409  
    410 !------------------------------------------------------------------------------!
     399
     400
     401!--------------------------------------------------------------------------------------------------!
    411402! Description:
    412403! ------------
    413 !> Calculation of the plant canopy transpiration rate based on the Jarvis-Stewart
    414 !> with parametrizations described in Daudet et al. (1999; Agricult. and Forest
    415 !> Meteorol. 97) and Ngao, Adam and Saudreau (2017;  Agricult. and Forest Meteorol
    416 !> 237-238). Model functions f1-f4 were adapted from Stewart (1998; Agric.
    417 !> and Forest. Meteorol. 43) instead, because they are valid for broader intervals
    418 !> of values. Funcion f4 used in form present in van Wijk et al. (1998;
    419 !> Tree Physiology 20).
     404!> Calculation of the plant canopy transpiration rate based on the Jarvis-Stewart with
     405!> parametrizations described in Daudet et al. (1999; Agricult. and Forest Meteorol. 97) and Ngao,
     406!> Adam and Saudreau (2017;  Agricult. and Forest Meteorol 237-238). Model functions f1-f4 were
     407!> adapted from Stewart (1998; Agric. and Forest. Meteorol. 43) instead, because they are valid for
     408!> broader intervals of values. Funcion f4 used in form present in van Wijk et al. (1998; Tree
     409!> Physiology 20).
    420410!>
    421 !> This subroutine is called from subroutine radiation_interaction
    422 !> after the calculation of radiation in plant canopy boxes.
     411!> This subroutine is called from subroutine radiation_interaction after the calculation of
     412!> radiation in plant canopy boxes.
    423413!> (arrays pcbinsw and pcbinlw).
    424414!>
    425 !------------------------------------------------------------------------------!
     415!--------------------------------------------------------------------------------------------------!
    426416 SUBROUTINE pcm_calc_transpiration_rate(i, j, k, kk, pcbsw, pcblw, pcbtr, pcblh)
    427417
    428418!
    429 !--  input parameters
     419!--  Input parameters
    430420     INTEGER(iwp), INTENT(IN) ::  i, j, k, kk        !< indices of the pc gridbox
     421     REAL(wp), INTENT(IN)     ::  pcblw              !< lw radiation in gridbox (W)
    431422     REAL(wp), INTENT(IN)     ::  pcbsw              !< sw radiation in gridbox (W)
    432      REAL(wp), INTENT(IN)     ::  pcblw              !< lw radiation in gridbox (W)
     423     REAL(wp), INTENT(OUT)    ::  pcblh              !< latent heat from transpiration dT/dt (K/s)
    433424     REAL(wp), INTENT(OUT)    ::  pcbtr              !< transpiration rate dq/dt (kg/kg/s)
    434      REAL(wp), INTENT(OUT)    ::  pcblh              !< latent heat from transpiration dT/dt (K/s)
    435 
    436 !--  variables and parameters for calculation of transpiration rate
    437      REAL(wp)             ::  sat_press, sat_press_d, temp, v_lad
    438      REAL(wp)             ::  d_fact, g_b, g_s, wind_speed, evapor_rate
    439      REAL(wp)             ::  f1, f2, f3, f4, vpd, rswc, e_eq, e_imp, rad
    440      REAL(wp), PARAMETER  ::  gama_psychr = 66.0_wp !< psychrometric constant (Pa/K)
    441      REAL(wp), PARAMETER  ::  g_s_max = 0.01        !< maximum stomatal conductivity (m/s)
    442      REAL(wp), PARAMETER  ::  m_soil = 0.4_wp       !< soil water content (needs to adjust or take from LSM)
    443      REAL(wp), PARAMETER  ::  m_wilt = 0.01_wp      !< wilting point soil water content (needs to adjust or take from LSM)
    444      REAL(wp), PARAMETER  ::  m_sat = 0.51_wp       !< saturation soil water content (needs to adjust or take from LSM)
    445      REAL(wp), PARAMETER  ::  t2_min = 0.0_wp       !< minimal temperature for calculation of f2
    446      REAL(wp), PARAMETER  ::  t2_max = 40.0_wp      !< maximal temperature for calculation of f2
    447 
    448 
     425
     426!--  Variables and parameters for calculation of transpiration rate
     427     REAL(wp), PARAMETER ::  gama_psychr = 66.0_wp !< psychrometric constant (Pa/K)
     428     REAL(wp), PARAMETER ::  g_s_max = 0.01        !< maximum stomatal conductivity (m/s)
     429     REAL(wp), PARAMETER ::  m_soil = 0.4_wp       !< soil water content (needs to adjust or take from LSM)
     430     REAL(wp), PARAMETER ::  m_wilt = 0.01_wp      !< wilting point soil water content (needs to adjust or take from LSM)
     431     REAL(wp), PARAMETER ::  m_sat = 0.51_wp       !< saturation soil water content (needs to adjust or take from LSM)
     432     REAL(wp), PARAMETER ::  t2_min = 0.0_wp       !< minimal temperature for calculation of f2
     433     REAL(wp), PARAMETER ::  t2_max = 40.0_wp      !< maximal temperature for calculation of f2
     434
     435     REAL(wp) ::  d_fact
     436     REAL(wp) ::  e_eq
     437     REAL(wp) ::  e_imp
     438     REAL(wp) ::  evapor_rate
     439     REAL(wp) ::  f1
     440     REAL(wp) ::  f2
     441     REAL(wp) ::  f3
     442     REAL(wp) ::  f4
     443     REAL(wp) ::  g_b
     444     REAL(wp) ::  g_s
     445     REAL(wp) ::  rad
     446     REAL(wp) ::  rswc
     447     REAL(wp) ::  sat_press
     448     REAL(wp) ::  sat_press_d
     449     REAL(wp) ::  temp
     450     REAL(wp) ::  v_lad
     451     REAL(wp) ::  vpd
     452     REAL(wp) ::  wind_speed
     453
     454!
    449455!--  Temperature (deg C)
    450456     temp = pt(k,j,i) * exner(k) - degc_to_k
     457!
    451458!--  Coefficient for conversion of radiation to grid to radiation to unit leaves surface
    452459     v_lad = 1.0_wp / ( MAX( lad_s(kk,j,i), 1.0E-10_wp ) * dx * dy * dz(1) )
     460!
    453461!--  Magnus formula for the saturation pressure (see Ngao, Adam and Saudreau (2017) eq. 1)
    454462!--  There are updated formulas available, kept consistent with the rest of the parametrization
    455      sat_press = 610.8_wp * exp(17.27_wp * temp/(temp + 237.3_wp))
     463     sat_press = 610.8_wp * EXP( 17.27_wp * temp / ( temp + 237.3_wp ) )
     464!
    456465!--  Saturation pressure derivative (derivative of the above)
    457      sat_press_d = sat_press * 17.27_wp * 237.3_wp / (temp + 237.3_wp)**2
     466     sat_press_d = sat_press * 17.27_wp * 237.3_wp / ( temp + 237.3_wp )**2
     467!
    458468!--  Wind speed
    459      wind_speed = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 +            &
    460                         ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 +            &
     469     wind_speed = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 +                                &
     470                        ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 +                                &
    461471                        ( 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) )**2 )
     472!
    462473!--  Aerodynamic conductivity (Daudet et al. (1999) eq. 14
    463474     g_b = 0.01_wp * wind_speed + 0.0071_wp
     475!
    464476!--  Radiation flux per leaf surface unit
    465477     rad = pcbsw * v_lad
     478!
    466479!--  First function for calculation of stomatal conductivity (radiation dependency)
    467480!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 17
    468      f1 = rad * (1000.0_wp+42.1_wp) / 1000.0_wp / (rad+42.1_wp)
     481     f1 = rad * ( 1000.0_wp + 42.1_wp ) / 1000.0_wp / ( rad + 42.1_wp )
     482!
    469483!--  Second function for calculation of stomatal conductivity (temperature dependency)
    470484!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 21
    471      f2 = MAX(t2_min, (temp-t2_min) * MAX(0.0_wp,t2_max-temp)**((t2_max-16.9_wp)/(16.9_wp-t2_min)) / &
    472               ((16.9_wp-t2_min) * (t2_max-16.9_wp)**((t2_max-16.9_wp)/(16.9_wp-t2_min))) )
     485     f2 = MAX( t2_min, ( temp - t2_min ) * MAX( 0.0_wp, t2_max - temp )**( ( t2_max - 16.9_wp ) /  &
     486                                                                           ( 16.9_wp - t2_min ) )  &
     487               / ( ( 16.9_wp - t2_min ) * ( t2_max - 16.9_wp )**( ( t2_max - 16.9_wp ) /           &
     488                                                                  ( 16.9_wp - t2_min ) ) ) )
     489!
    473490!--  Water pressure deficit
    474491!--  Ngao, Adam and Saudreau (2017) eq. 6 but with water vapour partial pressure
    475      vpd = max( sat_press - q(k,j,i) * hyp(k) / rd_d_rv, 0._wp )
     492     vpd = MAX( sat_press - q(k,j,i) * hyp(k) / rd_d_rv, 0._wp )
     493!
    476494!--  Third function for calculation of stomatal conductivity (water pressure deficit dependency)
    477495!--  Ngao, Adam and Saudreau (2017) Table 1, limited from below according to Stewart (1988)
    478 !--  The coefficients of the linear dependence should better correspond to broad-leaved trees
    479 !--  than the coefficients from Stewart (1988) which correspond to conifer trees.
    480      vpd = MIN(MAX(vpd,770.0_wp),3820.0_wp)
     496!--  The coefficients of the linear dependence should better correspond to broad-leaved trees than
     497!--  the coefficients from Stewart (1988) which correspond to conifer trees.
     498     vpd = MIN( MAX( vpd, 770.0_wp ), 3820.0_wp )
    481499     f3 = -2E-4_wp * vpd + 1.154_wp
     500!
    482501!--  Fourth function for calculation of stomatal conductivity (soil moisture dependency)
    483502!--  Residual soil water content
     
    485504!--  TODO - over LSM surface might be calculated from LSM parameters
    486505     rswc = ( m_sat - m_soil ) / ( m_sat - m_wilt )
    487 !--  van Wijk et al. (1998; Tree Physiology 20) eq. 5-6 (it is a reformulation of eq. 22-23 of Stewart(1988))
    488      f4 = MAX(0.0_wp, MIN(1.0_wp - 0.041_wp * EXP(3.2_wp * rswc), 1.0_wp - 0.041_wp))
     506!
     507!--  van Wijk et al. (1998; Tree Physiology 20) eq. 5-6 (it is a reformulation of eq. 22-23 of
     508!--  Stewart(1988))
     509     f4 = MAX( 0.0_wp, MIN( 1.0_wp - 0.041_wp * EXP( 3.2_wp * rswc ), 1.0_wp - 0.041_wp ) )
     510!
    489511!--  Stomatal conductivity
    490512!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 12
    491513!--  (notation according to Ngao, Adam and Saudreau (2017) and others)
    492514     g_s = g_s_max * f1 * f2 * f3 * f4 + 1.0E-10_wp
     515!
    493516!--  Decoupling factor
    494517!--  Daudet et al. (1999) eq. 6
    495      d_fact = (sat_press_d / gama_psychr + 2.0_wp ) /                        &
    496               (sat_press_d / gama_psychr + 2.0_wp + 2.0_wp * g_b / g_s )
     518     d_fact = ( sat_press_d / gama_psychr + 2.0_wp ) /                                             &
     519              ( sat_press_d / gama_psychr + 2.0_wp + 2.0_wp * g_b / g_s )
     520!
    497521!--  Equilibrium evaporation rate
    498522!--  Daudet et al. (1999) eq. 4
    499      e_eq = (pcbsw + pcblw) * v_lad * sat_press_d /                         &
    500                  gama_psychr /( sat_press_d / gama_psychr + 2.0_wp ) / l_v
     523     e_eq = ( pcbsw + pcblw ) * v_lad * sat_press_d /                                              &
     524            gama_psychr / ( sat_press_d / gama_psychr + 2.0_wp ) / l_v
     525!
    501526!--  Imposed evaporation rate
    502527!--  Daudet et al. (1999) eq. 5
    503528     e_imp = r_d * pt(k,j,i) * exner(k) / hyp(k) * c_p * g_s * vpd / gama_psychr / l_v
     529!
    504530!--  Evaporation rate
    505531!--  Daudet et al. (1999) eq. 3
    506532!--  (evaporation rate is limited to non-negative values)
    507      evapor_rate = MAX(d_fact * e_eq + ( 1.0_wp - d_fact ) * e_imp, 0.0_wp)
     533     evapor_rate = MAX( d_fact * e_eq + ( 1.0_wp - d_fact ) * e_imp, 0.0_wp )
     534!
    508535!--  Conversion of evaporation rate to q tendency in gridbox
    509536!--  dq/dt = E * LAD * V_g / (rho_air * V_g)
    510537     pcbtr = evapor_rate * r_d * pt(k,j,i) * exner(k) * lad_s(kk,j,i) / hyp(k)  !-- = dq/dt
     538!
    511539!--  latent heat from evaporation
    512540     pcblh = pcbtr * lv_d_cp  !-- = - dT/dt
     
    515543
    516544
    517 !------------------------------------------------------------------------------!
     545!--------------------------------------------------------------------------------------------------!
    518546! Description:
    519547! ------------
    520548!> Check data output for plant canopy model
    521 !------------------------------------------------------------------------------!
     549!--------------------------------------------------------------------------------------------------!
    522550 SUBROUTINE pcm_check_data_output( var, unit )
    523551
    524     CHARACTER (LEN=*) ::  unit  !< 
     552    CHARACTER (LEN=*) ::  unit  !<
    525553    CHARACTER (LEN=*) ::  var   !<
    526554
     
    530558       CASE ( 'pcm_heatrate' )
    531559!
    532 !--       Output of heatrate can be only done if it is explicitely set by cthf,
    533 !--       or parametrized by absorption of radiation. The latter, however, is
    534 !--       only available if radiation_interactions are on. Note, these are
    535 !--       enabled if land-surface or urban-surface is switched-on. Using
    536 !--       radiation_interactions_on directly is not possible since it belongs
    537 !--       to the radition_model, which in turn depends on the plant-canopy model,
    538 !--       creating circular dependencies.
    539           IF ( cthf == 0.0_wp  .AND. (  .NOT.  urban_surface  .AND.             &
    540                                         .NOT.  land_surface ) )  THEN
    541              message_string = 'output of "' // TRIM( var ) // '" requi' //      &
     560!--       Output of heatrate can be only done if it is explicitely set by cthf, or parametrized by
     561!--       absorption of radiation. The latter, however, is only available if radiation_interactions
     562!--       are on. Note, these are enabled if land-surface or urban-surface is switched-on. Using
     563!--       radiation_interactions_on directly is not possible since it belongs to the
     564!--       radition_model, which in turn depends on the plant-canopy model, creating circular
     565!--       dependencies.
     566          IF ( cthf == 0.0_wp  .AND.  ( .NOT. urban_surface .AND. .NOT. land_surface ) )  THEN
     567             message_string = 'output of "' // TRIM( var ) // '" requi' //                         &
    542568                              'res setting of parameter cthf /= 0.0'
    543569             CALL message( 'pcm_check_data_output', 'PA0718', 1, 2, 0, 6, 0 )
    544570          ENDIF
    545571          unit = 'K s-1'
    546    
     572
    547573       CASE ( 'pcm_transpirationrate' )
    548574          unit = 'kg kg-1 s-1'
     
    562588
    563589 END SUBROUTINE pcm_check_data_output
    564  
    565  
    566 !------------------------------------------------------------------------------!
     590
     591
     592!--------------------------------------------------------------------------------------------------!
    567593! Description:
    568594! ------------
    569595!> Check parameters routine for plant canopy model
    570 !------------------------------------------------------------------------------!
     596!--------------------------------------------------------------------------------------------------!
    571597    SUBROUTINE pcm_check_parameters
    572            
     598
    573599       IF ( ocean_mode )  THEN
    574           message_string = 'plant_canopy = .TRUE. is not allowed in the '//    &
    575                            'ocean'
     600          message_string = 'plant_canopy = .TRUE. is not allowed in the ocean'
    576601          CALL message( 'pcm_check_parameters', 'PA0696', 1, 2, 0, 6, 0 )
    577602       ENDIF
    578            
     603
    579604       IF ( canopy_drag_coeff == 0.0_wp )  THEN
    580           message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// &
     605          message_string = 'plant_canopy = .TRUE. requires a non-zero drag ' //                    &
    581606                           'coefficient & given value is canopy_drag_coeff = 0.0'
    582607          CALL message( 'pcm_check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
    583608       ENDIF
    584609
    585        IF ( ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad == 9999999.9_wp ) .OR.&
    586               beta_lad /= 9999999.9_wp   .AND.  alpha_lad == 9999999.9_wp )  THEN
    587           message_string = 'using the beta function for the construction ' //  &
    588                            'of the leaf area density profile requires '    //  &
     610       IF ( ( alpha_lad /= 9999999.9_wp .AND. beta_lad == 9999999.9_wp )  .OR.                     &
     611            beta_lad /= 9999999.9_wp  .AND.  alpha_lad == 9999999.9_wp )  THEN
     612          message_string = 'using the beta function for the construction ' //                      &
     613                           'of the leaf area density profile requires '    //                      &
    589614                           'both alpha_lad and beta_lad to be /= 9999999.9'
    590615          CALL message( 'pcm_check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
     
    592617
    593618       IF ( calc_beta_lad_profile  .AND.  lai_beta == 0.0_wp )  THEN
    594           message_string = 'using the beta function for the construction ' //  &
    595                            'of the leaf area density profile requires '    //  &
    596                            'a non-zero lai_beta, but given value is '      //  &
     619          message_string = 'using the beta function for the construction ' //                      &
     620                           'of the leaf area density profile requires '    //                      &
     621                           'a non-zero lai_beta, but given value is '      //                      &
    597622                           'lai_beta = 0.0'
    598623          CALL message( 'pcm_check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
     
    600625
    601626       IF ( calc_beta_lad_profile  .AND.  lad_surface /= 0.0_wp )  THEN
    602           message_string = 'simultaneous setting of alpha_lad /= 9999999.9 '// &
    603                            'combined with beta_lad /= 9999999.9 '           // &
    604                            'and lad_surface /= 0.0 is not possible, '       // &
    605                            'use either vertical gradients or the beta '     // &
    606                            'function for the construction of the leaf area '// &
     627          message_string = 'simultaneous setting of alpha_lad /= 9999999.9 '//                     &
     628                           'combined with beta_lad /= 9999999.9 '           //                     &
     629                           'and lad_surface /= 0.0 is not possible, '       //                     &
     630                           'use either vertical gradients or the beta '     //                     &
     631                           'function for the construction of the leaf area '//                     &
    607632                           'density profile'
    608633          CALL message( 'pcm_check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
    609        ENDIF 
     634       ENDIF
    610635
    611636       IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
    612           message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
    613                           ' seifert_beheng'
     637          message_string = 'plant_canopy = .TRUE. requires cloud_scheme /= seifert_beheng'
    614638          CALL message( 'pcm_check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
    615639       ENDIF
    616640
    617     END SUBROUTINE pcm_check_parameters 
    618  
    619 
    620 !------------------------------------------------------------------------------!
     641    END SUBROUTINE pcm_check_parameters
     642
     643
     644!--------------------------------------------------------------------------------------------------!
    621645!
    622646! Description:
    623647! ------------
    624648!> Subroutine for averaging 3D data
    625 !------------------------------------------------------------------------------!
     649!--------------------------------------------------------------------------------------------------!
    626650 SUBROUTINE pcm_3d_data_averaging( mode, variable )
    627651
    628652    CHARACTER (LEN=*) ::  mode    !<
    629     CHARACTER (LEN=*) :: variable !<
     653    CHARACTER (LEN=*) ::  variable !<
    630654
    631655    INTEGER(iwp) ::  i            !<
     
    663687       END SELECT
    664688
    665     ELSEIF ( mode == 'sum' )  THEN
     689    ELSE IF ( mode == 'sum' )  THEN
    666690
    667691       SELECT CASE ( TRIM( variable ) )
    668692
    669693          CASE ( 'pcm_heatrate' )
    670              IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
     694             IF ( ALLOCATED( pcm_heatrate_av ) )  THEN
    671695                DO  i = nxl, nxr
    672696                   DO  j = nys, nyn
    673697                      IF ( pch_index_ji(j,i) /= 0 )  THEN
    674698                         DO  k = 0, pch_index_ji(j,i)
    675                             pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i) + pcm_heating_rate(k,j,i)
     699                            pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i) +                      &
     700                                                     pcm_heating_rate(k,j,i)
    676701                         ENDDO
    677702                      ENDIF
     
    682707
    683708          CASE ( 'pcm_latentrate' )
    684              IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
     709             IF ( ALLOCATED( pcm_latentrate_av ) )  THEN
    685710                DO  i = nxl, nxr
    686711                   DO  j = nys, nyn
    687712                      IF ( pch_index_ji(j,i) /= 0 )  THEN
    688713                         DO  k = 0, pch_index_ji(j,i)
    689                             pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i) + pcm_latent_rate(k,j,i)
     714                            pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i) +                  &
     715                                                       pcm_latent_rate(k,j,i)
    690716                         ENDDO
    691717                      ENDIF
     
    696722
    697723          CASE ( 'pcm_transpirationrate' )
    698              IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
     724             IF ( ALLOCATED( pcm_transpirationrate_av ) )  THEN
    699725                DO  i = nxl, nxr
    700726                   DO  j = nys, nyn
    701727                      IF ( pch_index_ji(j,i) /= 0 )  THEN
    702728                         DO  k = 0, pch_index_ji(j,i)
    703                             pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i) + pcm_transpiration_rate(k,j,i)
     729                            pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i) +    &
     730                                                              pcm_transpiration_rate(k,j,i)
    704731                         ENDDO
    705732                      ENDIF
     
    713740       END SELECT
    714741
    715     ELSEIF ( mode == 'average' )  THEN
     742    ELSE IF ( mode == 'average' )  THEN
    716743
    717744       SELECT CASE ( TRIM( variable ) )
    718745
    719746          CASE ( 'pcm_heatrate' )
    720              IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
     747             IF ( ALLOCATED( pcm_heatrate_av ) )  THEN
    721748                DO  i = nxlg, nxrg
    722749                   DO  j = nysg, nyng
    723750                      IF ( pch_index_ji(j,i) /= 0 )  THEN
    724751                         DO  k = 0, pch_index_ji(j,i)
    725                             pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i)                 &
     752                            pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i)                        &
    726753                                                     / REAL( average_count_3d, KIND=wp )
    727754                         ENDDO
     
    733760
    734761          CASE ( 'pcm_latentrate' )
    735              IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
     762             IF ( ALLOCATED( pcm_latentrate_av ) )  THEN
    736763                DO  i = nxlg, nxrg
    737764                   DO  j = nysg, nyng
    738765                      IF ( pch_index_ji(j,i) /= 0 )  THEN
    739766                         DO  k = 0, pch_index_ji(j,i)
    740                             pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i)              &
     767                            pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i)                    &
    741768                                                       / REAL( average_count_3d, KIND=wp )
    742769                         ENDDO
     
    748775
    749776          CASE ( 'pcm_transpirationrate' )
    750              IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
     777             IF ( ALLOCATED( pcm_transpirationrate_av ) )  THEN
    751778                DO  i = nxlg, nxrg
    752779                   DO  j = nysg, nyng
    753780                      IF ( pch_index_ji(j,i) /= 0 )  THEN
    754781                         DO  k = 0, pch_index_ji(j,i)
    755                             pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i)  &
     782                            pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i)      &
    756783                                                              / REAL( average_count_3d, KIND=wp )
    757784                         ENDDO
     
    767794 END SUBROUTINE pcm_3d_data_averaging
    768795
    769 !------------------------------------------------------------------------------!
     796!--------------------------------------------------------------------------------------------------!
    770797!
    771798! Description:
    772799! ------------
    773 !> Subroutine defining 3D output variables.
    774 !> Note, 3D plant-canopy output has it's own vertical output dimension, meaning
    775 !> that 3D output is relative to the model surface now rather than at the actual
    776 !> grid point where the plant canopy is located.
    777 !------------------------------------------------------------------------------!
    778  SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf, fill_value,     &
    779                                 nzb_do, nzt_do )
     800!> Subroutine defining 3D output variables.
     801!> Note, 3D plant-canopy output has it's own vertical output dimension, meaning that 3D output is
     802!> relative to the model surface now rather than at the actual grid point where the plant canopy is
     803!> located.
     804!--------------------------------------------------------------------------------------------------!
     805 SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
    780806
    781807    CHARACTER (LEN=*) ::  variable !< treated variable
     
    791817
    792818    REAL(wp)     ::  fill_value !< fill value
     819
    793820    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< data output array
    794821
     
    801828!
    802829!--    Note, to save memory arrays for heating are allocated from 0:pch_index.
    803 !--    Thus, output must be relative to these array indices. Further, check
    804 !--    whether the output is within the vertical output range,
    805 !--    i.e. nzb_do:nzt_do, which is necessary as local_pf is only allocated
    806 !--    for this index space. Note, plant-canopy output has a separate
    807 !--    vertical output coordinate zlad, so that output is mapped down to the
    808 !--    surface.
     830!--    Thus, output must be relative to these array indices. Further, check whether the output is
     831!--    within the vertical output range, i.e. nzb_do:nzt_do, which is necessary as local_pf is only
     832!--    allocated for this index space. Note, plant-canopy output has a separate vertical output
     833!--    coordinate zlad, so that output is mapped down to the surface.
    809834       CASE ( 'pcm_heatrate' )
    810835          IF ( av == 0 )  THEN
     
    891916    END SELECT
    892917
    893  END SUBROUTINE pcm_data_output_3d 
    894          
    895 !------------------------------------------------------------------------------!
     918 END SUBROUTINE pcm_data_output_3d
     919
     920!--------------------------------------------------------------------------------------------------!
    896921!
    897922! Description:
     
    899924!> Subroutine defining appropriate grid for netcdf variables.
    900925!> It is called from subroutine netcdf.
    901 !------------------------------------------------------------------------------!
     926!--------------------------------------------------------------------------------------------------!
    902927 SUBROUTINE pcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
    903928
    904      CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
    905      LOGICAL, INTENT(OUT)           ::  found       !<
    906      CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
    907      CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
    908      CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
     929     CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
     930     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
     931     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
     932     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
     933
     934     LOGICAL, INTENT(OUT)           ::  found       !<
     935
    909936
    910937     found  = .TRUE.
     
    914941     SELECT CASE ( TRIM( var ) )
    915942
    916         CASE ( 'pcm_heatrate', 'pcm_bad', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate')
     943        CASE ( 'pcm_heatrate', 'pcm_bad', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate' )
    917944           grid_x = 'x'
    918945           grid_y = 'y'
     
    927954
    928955 END SUBROUTINE pcm_define_netcdf_grid
    929  
    930  
    931 !------------------------------------------------------------------------------!
     956
     957
     958!--------------------------------------------------------------------------------------------------!
    932959! Description:
    933960! ------------
    934961!> Header output for plant canopy model
    935 !------------------------------------------------------------------------------!
     962!--------------------------------------------------------------------------------------------------!
    936963 SUBROUTINE pcm_header ( io )
    937  
     964
    938965    CHARACTER (LEN=10) ::  coor_chr            !<
    939 
    940966    CHARACTER (LEN=86) ::  coordinates         !<
    941967    CHARACTER (LEN=86) ::  gradients           !<
    942968    CHARACTER (LEN=86) ::  leaf_area_density   !<
    943969    CHARACTER (LEN=86) ::  slices              !<
    944  
    945     INTEGER(iwp) :: i                !<
    946     INTEGER(iwp),  INTENT(IN) ::  io !< Unit of the output file
    947     INTEGER(iwp) :: k                !<       
    948  
     970
     971    INTEGER(iwp)              ::  i   !<
     972    INTEGER(iwp),  INTENT(IN) ::  io  !< Unit of the output file
     973    INTEGER(iwp)              ::  k   !<
     974
    949975    REAL(wp) ::  canopy_height       !< canopy height (in m)
    950    
     976
     977
    951978    canopy_height = zw(pch_index)
    952979
    953     WRITE ( io, 1 )  canopy_mode, canopy_height, pch_index,                    &
    954                        canopy_drag_coeff                                       
    955     IF ( passive_scalar )  THEN                                               
    956        WRITE ( io, 2 )  leaf_scalar_exch_coeff,                                &
    957                           leaf_surface_conc
     980    WRITE( io, 1 )  canopy_mode, canopy_height, pch_index, canopy_drag_coeff
     981    IF ( passive_scalar )  THEN
     982       WRITE( io, 2 )  leaf_scalar_exch_coeff, leaf_surface_conc
    958983    ENDIF
    959984
    960985!
    961 !  Heat flux at the top of vegetation
    962     WRITE ( io, 3 )  cthf
    963 
    964 !
    965 !   Leaf area density profile, calculated either from given vertical
    966 !   gradients or from beta probability density function.
    967     IF (  .NOT. calc_beta_lad_profile )  THEN
     986!-- Heat flux at the top of vegetation
     987    WRITE( io, 3 )  cthf
     988
     989!
     990!-- Leaf area density profile, calculated either from given vertical gradients or from beta
     991!-- probability density function.
     992    IF ( .NOT. calc_beta_lad_profile )  THEN
    968993
    969994!      Building output strings, starting with surface value
    970        WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
     995       WRITE( leaf_area_density, '(F7.4)' )  lad_surface
    971996       gradients = '------'
    972997       slices = '     0'
    973998       coordinates = '   0.0'
    974        DO i = 1, UBOUND(lad_vertical_gradient_level_ind, DIM=1)
    975           IF  ( lad_vertical_gradient_level_ind(i) /= -9999 ) THEN
    976 
    977              WRITE (coor_chr,'(F7.2)') lad(lad_vertical_gradient_level_ind(i))
     999       DO  i = 1, UBOUND( lad_vertical_gradient_level_ind, DIM=1 )
     1000          IF  ( lad_vertical_gradient_level_ind(i) /= -9999 )  THEN
     1001
     1002             WRITE( coor_chr, '(F7.2)' ) lad(lad_vertical_gradient_level_ind(i))
    9781003             leaf_area_density = TRIM( leaf_area_density ) // ' ' // TRIM( coor_chr )
    9791004
    980              WRITE (coor_chr,'(F7.2)') lad_vertical_gradient(i)
     1005             WRITE( coor_chr, '(F7.2)' ) lad_vertical_gradient(i)
    9811006             gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
    9821007
    983              WRITE (coor_chr,'(I7)') lad_vertical_gradient_level_ind(i)
     1008             WRITE( coor_chr, '(I7)' ) lad_vertical_gradient_level_ind(i)
    9841009             slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
    9851010
    986              WRITE (coor_chr,'(F7.1)') lad_vertical_gradient_level(i)
     1011             WRITE( coor_chr, '(F7.1)' ) lad_vertical_gradient_level(i)
    9871012             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
    9881013          ELSE
     
    9911016       ENDDO
    9921017
    993        WRITE ( io, 4 )  TRIM( coordinates ), TRIM( leaf_area_density ),        &
    994                           TRIM( gradients ), TRIM( slices )
     1018       WRITE( io, 4 )  TRIM( coordinates ), TRIM( leaf_area_density ), TRIM( gradients ),          &
     1019                       TRIM( slices )
    9951020
    9961021    ELSE
    997    
    998        WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
     1022
     1023       WRITE( leaf_area_density, '(F7.4)' )  lad_surface
    9991024       coordinates = '   0.0'
    1000        
     1025
    10011026       DO  k = 1, pch_index
    10021027
    1003           WRITE (coor_chr,'(F7.2)')  lad(k)
    1004           leaf_area_density = TRIM( leaf_area_density ) // ' ' //              &
    1005                               TRIM( coor_chr )                                 
    1006                                                                                
    1007           WRITE (coor_chr,'(F7.1)')  zu(k)                                     
    1008           coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )       
    1009                                                                                
    1010        ENDDO                                                                   
    1011                                                                                
    1012        WRITE ( io, 5 ) TRIM( coordinates ), TRIM( leaf_area_density ),         &
    1013                        alpha_lad, beta_lad, lai_beta
    1014 
    1015     ENDIF 
    1016 
    1017 1   FORMAT (//' Vegetation canopy (drag) model:'/                              &
    1018               ' ------------------------------'//                              &
    1019               ' Canopy mode: ', A /                                            &
    1020               ' Canopy height: ',F6.2,'m (',I4,' grid points)' /               &
    1021               ' Leaf drag coefficient: ',F6.2 /)
    1022 2   FORMAT (/ ' Scalar exchange coefficient: ',F6.2 /                          &
    1023               ' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /)
    1024 3   FORMAT (' Predefined constant heatflux at the top of the vegetation: ',    &
    1025              F6.2, ' K m/s')
    1026 4   FORMAT (/ ' Characteristic levels of the leaf area density:'//             &
    1027               ' Height:              ',A,'  m'/                                &
    1028               ' Leaf area density:   ',A,'  m**2/m**3'/                        &
    1029               ' Gradient:            ',A,'  m**2/m**4'/                        &
    1030               ' Gridpoint:           ',A)
    1031 5   FORMAT (//' Characteristic levels of the leaf area density and coefficients:'&
    1032           //  ' Height:              ',A,'  m'/                                &
    1033               ' Leaf area density:   ',A,'  m**2/m**3'/                        &
    1034               ' Coefficient alpha: ',F6.2 /                                    &
    1035               ' Coefficient beta: ',F6.2 /                                     &
    1036               ' Leaf area index: ',F6.2,'  m**2/m**2' /)   
    1037        
     1028          WRITE( coor_chr,'(F7.2)' )  lad(k)
     1029          leaf_area_density = TRIM( leaf_area_density ) // ' ' // TRIM( coor_chr )
     1030
     1031          WRITE(coor_chr,'(F7.1)')  zu(k)
     1032          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
     1033
     1034       ENDDO
     1035
     1036       WRITE( io, 5 ) TRIM( coordinates ), TRIM( leaf_area_density ), alpha_lad, beta_lad, lai_beta
     1037
     1038    ENDIF
     1039
     10401   FORMAT (/ /' Vegetation canopy (drag) model:' / ' ------------------------------' //           &
     1041              ' Canopy mode: ', A / ' Canopy height: ', F6.2, 'm (',I4,' grid points)' /           &
     1042              ' Leaf drag coefficient: ', F6.2 /)
     10432   FORMAT (/ ' Scalar exchange coefficient: ',F6.2 /                                              &
     1044              ' Scalar concentration at leaf surfaces in kg/m**3: ', F6.2 /)
     10453   FORMAT ( ' Predefined constant heatflux at the top of the vegetation: ', F6.2, ' K m/s')
     10464   FORMAT (/ ' Characteristic levels of the leaf area density:' //                                &
     1047              ' Height:              ', A, '  m' /                                                 &
     1048              ' Leaf area density:   ', A, '  m**2/m**3' /                                         &
     1049              ' Gradient:            ', A, '  m**2/m**4' /                                         &
     1050              ' Gridpoint:           ', A )
     10515   FORMAT (//' Characteristic levels of the leaf area density and coefficients:' //               &
     1052              ' Height:              ', A, '  m' /                                                 &
     1053              ' Leaf area density:   ', A, '  m**2/m**3' /                                         &
     1054              ' Coefficient alpha: ',F6.2 /                                                        &
     1055              ' Coefficient beta: ',F6.2 /                                                         &
     1056              ' Leaf area index: ',F6.2,'  m**2/m**2' /)
     1057
    10381058    END SUBROUTINE pcm_header
    1039  
    1040  
    1041 !------------------------------------------------------------------------------!
     1059
     1060
     1061!--------------------------------------------------------------------------------------------------!
    10421062! Description:
    10431063! ------------
    10441064!> Initialization of the plant canopy model
    1045 !------------------------------------------------------------------------------!
     1065!--------------------------------------------------------------------------------------------------!
    10461066    SUBROUTINE pcm_init
    10471067
    1048        USE exchange_horiz_mod,                                                    &
     1068       USE exchange_horiz_mod,                                                                     &
    10491069           ONLY:  exchange_horiz
    10501070
     
    10541074       INTEGER(iwp) ::  m   !< running index
    10551075
    1056        LOGICAL      ::  lad_on_top = .FALSE.  !< dummy flag to indicate that LAD is defined on a building roof
    1057        LOGICAL      ::  bad_on_top = .FALSE.  !< dummy flag to indicate that BAD is defined on a building roof
     1076       LOGICAL ::  lad_on_top = .FALSE.  !< dummy flag to indicate that LAD is defined on a building roof
     1077       LOGICAL ::  bad_on_top = .FALSE.  !< dummy flag to indicate that BAD is defined on a building roof
    10581078
    10591079       REAL(wp) ::  canopy_height   !< canopy height for lad-profile construction
    10601080       REAL(wp) ::  gradient        !< gradient for lad-profile construction
    1061        REAL(wp) ::  int_bpdf        !< vertical integral for lad-profile construction     
     1081       REAL(wp) ::  int_bpdf        !< vertical integral for lad-profile construction
    10621082       REAL(wp) ::  lad_max         !< maximum LAD value in the model domain, used to perform a check
    10631083
     1084
    10641085       IF ( debug_output )  CALL debug_message( 'pcm_init', 'start' )
    10651086!
    1066 !--    Allocate one-dimensional arrays for the computation of the
    1067 !--    leaf area density (lad) profile
     1087!--    Allocate one-dimensional arrays for the computation of the leaf area density (lad) profile
    10681088       ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
    10691089       lad = 0.0_wp
     
    10711091
    10721092!
    1073 !--    Set flag that indicates that the lad-profile shall be calculated by using
    1074 !--    a beta probability density function
     1093!--    Set flag that indicates that the lad-profile shall be calculated by using a beta probability
     1094!--    density function
    10751095       IF ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad /= 9999999.9_wp )  THEN
    10761096          calc_beta_lad_profile = .TRUE.
    10771097       ENDIF
    1078        
    1079        
    1080 !
    1081 !--    Compute the profile of leaf area density used in the plant
    1082 !--    canopy model. The profile can either be constructed from
    1083 !--    prescribed vertical gradients of the leaf area density or by
    1084 !--    using a beta probability density function (see e.g. Markkanen et al.,
    1085 !--    2003: Boundary-Layer Meteorology, 106, 437-459)
    1086        IF (  .NOT.  calc_beta_lad_profile )  THEN   
    1087 
    1088 !
    1089 !--       Use vertical gradients for lad-profile construction   
     1098
     1099
     1100!
     1101!--    Compute the profile of leaf area density used in the plant canopy model. The profile can
     1102!--    either be constructed from prescribed vertical gradients of the leaf area density or by using
     1103!--    a beta probability density function (see e.g. Markkanen et al., 2003: Boundary-Layer
     1104!--    Meteorology, 106, 437-459)
     1105       IF ( .NOT. calc_beta_lad_profile )  THEN
     1106
     1107!
     1108!--       Use vertical gradients for lad-profile construction
    10901109          i = 1
    10911110          gradient = 0.0_wp
     
    10941113          lad_vertical_gradient_level_ind(1) = 0
    10951114
    1096           DO k = 1, pch_index
     1115          DO  k = 1, pch_index
    10971116             IF ( i < 11 )  THEN
    1098                 IF ( lad_vertical_gradient_level(i) < zu(k)  .AND.          &
     1117                IF ( lad_vertical_gradient_level(i) < zu(k)  .AND.                                 &
    10991118                     lad_vertical_gradient_level(i) >= 0.0_wp )  THEN
    11001119                   gradient = lad_vertical_gradient(i)
     
    11151134
    11161135!
    1117 !--       In case of no given leaf area density gradients, choose a vanishing
    1118 !--       gradient. This information is used for the HEADER and the RUN_CONTROL
    1119 !--       file.
     1136!--       In case of no given leaf area density gradients, choose a vanishing gradient. This
     1137!--       information is used for the HEADER and the RUN_CONTROL file.
    11201138          IF ( lad_vertical_gradient_level(1) == -9999999.9_wp )  THEN
    11211139             lad_vertical_gradient_level(1) = 0.0_wp
     
    11241142       ELSE
    11251143
    1126 ! 
     1144!
    11271145!--       Use beta function for lad-profile construction
    11281146          int_bpdf = 0.0_wp
    11291147          canopy_height = zw(pch_index)
    11301148
    1131           DO k = 0, pch_index
    1132              int_bpdf = int_bpdf +                                             &
    1133                       ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) *  &
    1134                       ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(              &
    1135                           beta_lad-1.0_wp ) )                                  &
    1136                       * ( ( zw(k+1)-zw(k) ) / canopy_height ) )
     1149          DO  k = 0, pch_index
     1150             int_bpdf = int_bpdf +                                                                 &
     1151                                 ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) *           &
     1152                                   ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) &
     1153                                   * ( ( zw(k+1)-zw(k) ) / canopy_height ) )
    11371154          ENDDO
    11381155
    11391156!
    11401157!--       Preliminary lad profile (defined on w-grid)
    1141           DO k = 0, pch_index
    1142              pre_lad(k) =  lai_beta *                                          &
    1143                         ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) )  &
    1144                         * ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(          &
    1145                               beta_lad-1.0_wp ) ) / int_bpdf                   &
    1146                         ) / canopy_height
     1158          DO  k = 0, pch_index
     1159             pre_lad(k) =  lai_beta *                                                              &
     1160                              ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) )                &
     1161                                * ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) )  &
     1162                                / int_bpdf                                                         &
     1163                              ) / canopy_height
    11471164          ENDDO
    11481165
    11491166!
    1150 !--       Final lad profile (defined on scalar-grid level, since most prognostic
    1151 !--       quantities are defined there, hence, less interpolation is required
    1152 !--       when calculating the canopy tendencies)
     1167!--       Final lad profile (defined on scalar-grid level, since most prognostic quantities are
     1168!--       defined there, hence, less interpolation is required when calculating the canopy
     1169!--       tendencies)
    11531170          lad(0) = pre_lad(0)
    1154           DO k = 1, pch_index
     1171          DO  k = 1, pch_index
    11551172             lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
    11561173          ENDDO
     
    11591176
    11601177!
    1161 !--    Allocate 3D-array for the leaf-area density (lad_s) as well as
    1162 !--    for basal-area densitiy (bad_s). Note, by default bad_s is zero.
     1178!--    Allocate 3D-array for the leaf-area density (lad_s) as well as for basal-area densitiy
     1179!--    (bad_s). Note, by default bad_s is zero.
    11631180       ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    11641181       ALLOCATE( bad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     
    11661183!
    11671184!--    Initialization of the canopy coverage in the model domain:
    1168 !--    Setting the parameter canopy_mode = 'homogeneous' initializes a canopy, which
    1169 !--    fully covers the domain surface
     1185!--    Setting the parameter canopy_mode = 'homogeneous' initializes a canopy, which fully covers
     1186!--    the domain surface
    11701187       SELECT CASE ( TRIM( canopy_mode ) )
    11711188
    1172           CASE( 'homogeneous' )
     1189          CASE ( 'homogeneous' )
    11731190
    11741191             DO  i = nxlg, nxrg
     
    11851202!--             Open the static input file
    11861203#if defined( __netcdf )
    1187                 CALL open_read_file( TRIM( input_file_static ) //              &
    1188                                      TRIM( coupling_char ),                    &
     1204                CALL open_read_file( TRIM( input_file_static ) //                                  &
     1205                                     TRIM( coupling_char ),                                        &
    11891206                                     pids_id )
    11901207
     
    11981215                IF ( check_existence( vars_pids, 'lad' ) )  THEN
    11991216                   leaf_area_density_f%from_file = .TRUE.
    1200                    CALL get_attribute( pids_id, char_fill,                     &
    1201                                        leaf_area_density_f%fill,               &
     1217                   CALL get_attribute( pids_id, char_fill,                                         &
     1218                                       leaf_area_density_f%fill,                                   &
    12021219                                       .FALSE., 'lad' )
    12031220!
    12041221!--                Inquire number of vertical vegetation layer
    1205                    CALL get_dimension_length( pids_id,                         &
    1206                                               leaf_area_density_f%nz,          &
     1222                   CALL get_dimension_length( pids_id,                                             &
     1223                                              leaf_area_density_f%nz,                              &
    12071224                                              'zlad' )
    12081225!
    12091226!--                Allocate variable for leaf-area density
    1210                    ALLOCATE( leaf_area_density_f%var                           &
    1211                                                 (0:leaf_area_density_f%nz-1,   &
    1212                                                  nys:nyn,nxl:nxr) )
    1213 
    1214                    CALL get_variable( pids_id, 'lad', leaf_area_density_f%var, &
    1215                                       nxl, nxr, nys, nyn,                      &
     1227                   ALLOCATE( leaf_area_density_f%var                                               &
     1228                                                 (0:leaf_area_density_f%nz-1,nys:nyn,nxl:nxr) )
     1229
     1230                   CALL get_variable( pids_id, 'lad', leaf_area_density_f%var, nxl, nxr, nys, nyn, &
    12161231                                      0, leaf_area_density_f%nz-1 )
    12171232
     
    12231238                IF ( check_existence( vars_pids, 'bad' ) )  THEN
    12241239                   basal_area_density_f%from_file = .TRUE.
    1225                    CALL get_attribute( pids_id, char_fill,                     &
    1226                                        basal_area_density_f%fill,              &
     1240                   CALL get_attribute( pids_id, char_fill,                                         &
     1241                                       basal_area_density_f%fill,                                  &
    12271242                                       .FALSE., 'bad' )
    12281243!
    12291244!--                Inquire number of vertical vegetation layer
    1230                    CALL get_dimension_length( pids_id,                         &
    1231                                               basal_area_density_f%nz,         &
     1245                   CALL get_dimension_length( pids_id,                                             &
     1246                                              basal_area_density_f%nz,                             &
    12321247                                              'zlad' )
    12331248!
    12341249!--                Allocate variable
    1235                    ALLOCATE( basal_area_density_f%var                          &
    1236                                               (0:basal_area_density_f%nz-1,    &
    1237                                                nys:nyn,nxl:nxr) )
    1238 
    1239                    CALL get_variable( pids_id, 'bad', basal_area_density_f%var,&
    1240                                       nxl, nxr, nys, nyn,                      &
    1241                                       0,  basal_area_density_f%nz-1 )
     1250                   ALLOCATE( basal_area_density_f%var                                              &
     1251                                                  (0:basal_area_density_f%nz-1,nys:nyn,nxl:nxr) )
     1252
     1253                   CALL get_variable( pids_id, 'bad', basal_area_density_f%var, nxl, nxr, nys, nyn,&
     1254                                      0, basal_area_density_f%nz-1 )
    12421255                ELSE
    12431256                   basal_area_density_f%from_file = .FALSE.
     
    12471260                IF ( check_existence( vars_pids, 'root_area_dens_r' ) )  THEN
    12481261                   root_area_density_lad_f%from_file = .TRUE.
    1249                    CALL get_attribute( pids_id, char_fill,                     &
    1250                                        root_area_density_lad_f%fill,           &
     1262                   CALL get_attribute( pids_id, char_fill,                                         &
     1263                                       root_area_density_lad_f%fill,                               &
    12511264                                       .FALSE., 'root_area_dens_r' )
    12521265!
    12531266!--                Inquire number of vertical soil layers
    1254                    CALL get_dimension_length( pids_id,                         &
    1255                                               root_area_density_lad_f%nz,      &
     1267                   CALL get_dimension_length( pids_id,                                             &
     1268                                              root_area_density_lad_f%nz,                          &
    12561269                                              'zsoil' )
    12571270!
    12581271!--                Allocate variable
    1259                    ALLOCATE( root_area_density_lad_f%var                       &
    1260                                                (0:root_area_density_lad_f%nz-1,&
    1261                                                 nys:nyn,nxl:nxr) )
    1262 
    1263                    CALL get_variable( pids_id, 'root_area_dens_r',             &
    1264                                       root_area_density_lad_f%var,             &
    1265                                       nxl, nxr, nys, nyn,                      &
    1266                                       0,  root_area_density_lad_f%nz-1 )
     1272                   ALLOCATE( root_area_density_lad_f%var                                           &
     1273                                                 (0:root_area_density_lad_f%nz-1,nys:nyn,nxl:nxr) )
     1274
     1275                   CALL get_variable( pids_id, 'root_area_dens_r', root_area_density_lad_f%var,    &
     1276                                      nxl, nxr, nys, nyn, 0, root_area_density_lad_f%nz-1 )
    12671277                ELSE
    12681278                   root_area_density_lad_f%from_file = .FALSE.
     
    12771287
    12781288!
    1279 !--          Initialize LAD with data from file. If LAD is given in NetCDF file,
    1280 !--          use these values, else take LAD profiles from ASCII file.
    1281 !--          Please note, in NetCDF file LAD is only given up to the maximum
    1282 !--          canopy top, indicated by leaf_area_density_f%nz. 
     1289!--          Initialize LAD with data from file. If LAD is given in NetCDF file, use these values,
     1290!--          else take LAD profiles from ASCII file.
     1291!--          Please note, in NetCDF file LAD is only given up to the maximum canopy top, indicated
     1292!--          by leaf_area_density_f%nz.
    12831293             lad_s = 0.0_wp
    12841294             IF ( leaf_area_density_f%from_file )  THEN
    12851295!
    1286 !--             Set also pch_index, used to be the upper bound of the vertical
    1287 !--             loops. Therefore, use the global top of the canopy layer.
     1296!--             Set also pch_index, used to be the upper bound of the vertical loops. Therefore, use
     1297!--             the global top of the canopy layer.
    12881298                pch_index = leaf_area_density_f%nz - 1
    12891299
     
    12911301                   DO  j = nys, nyn
    12921302                      DO  k = 0, leaf_area_density_f%nz - 1
    1293                          IF ( leaf_area_density_f%var(k,j,i) /=                &
    1294                               leaf_area_density_f%fill )                       &
    1295                          lad_s(k,j,i) = leaf_area_density_f%var(k,j,i)
     1303                         IF ( leaf_area_density_f%var(k,j,i) /=  leaf_area_density_f%fill )        &
     1304                            lad_s(k,j,i) = leaf_area_density_f%var(k,j,i)
    12961305                      ENDDO
    12971306!
    12981307!--                   Check if resolved vegetation is mapped onto buildings.
    1299 !--                   In general, this is allowed and also meaningful, e.g.
    1300 !--                   when trees carry across roofs. However, due to the
    1301 !--                   topography filtering, new building grid points can emerge
    1302 !--                   at location where also plant canopy is defined. As a
    1303 !--                   result, plant canopy is mapped on top of roofs, with
    1304 !--                   siginficant impact on the downstream flow field and the
    1305 !--                   nearby surface radiation. In order to avoid that
    1306 !--                   plant canopy is mistakenly mapped onto building roofs,
    1307 !--                   check for building grid points (bit 6) that emerge from
    1308 !--                   the filtering (bit 4) and set LAD to zero at these
    1309 !--                   artificially  created building grid points. This case,
    1310 !--                   an informative message is given.
    1311                       IF ( ANY( lad_s(:,j,i) /= 0.0_wp )  .AND.                &
    1312                            ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND.  &
     1308!--                   In general, this is allowed and also meaningful, e.g. when trees carry across
     1309!--                   roofs. However, due to the topography filtering, new building grid points can
     1310!--                   emerge at locations where also plant canopy is defined. As a result, plant
     1311!--                   canopy is mapped on top of roofs, with siginficant impact on the downstream
     1312!--                   flow field and the nearby surface radiation. In order to avoid that plant
     1313!--                   canopy is mistakenly mapped onto building roofs, check for building grid
     1314!--                   points (bit 6) that emerge from the filtering (bit 4) and set LAD to zero at
     1315!--                   these artificially created building grid points. This case, an informative
     1316!--                   message is given.
     1317                      IF ( ANY( lad_s(:,j,i) /= 0.0_wp )                 .AND.                     &
     1318                           ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) )  .AND.                     &
    13131319                           ANY( BTEST( wall_flags_total_0(:,j,i), 4 ) ) )  THEN
    13141320                         lad_s(:,j,i) = 0.0_wp
     
    13181324                ENDDO
    13191325#if defined( __parallel )
    1320                CALL MPI_ALLREDUCE( MPI_IN_PLACE, lad_on_top, 1, MPI_LOGICAL,  &
    1321                                    MPI_LOR, comm2d, ierr)
     1326               CALL MPI_ALLREDUCE( MPI_IN_PLACE, lad_on_top, 1, MPI_LOGICAL, MPI_LOR, comm2d, ierr)
    13221327#endif
    13231328                IF ( lad_on_top )  THEN
     
    13321337!
    13331338!            ASCII file
    1334 !--          Initialize canopy parameters canopy_drag_coeff,
    1335 !--          leaf_scalar_exch_coeff, leaf_surface_conc
    1336 !--          from file which contains complete 3D data (separate vertical profiles for
    1337 !--          each location).
     1339!--          Initialize canopy parameters canopy_drag_coeff, leaf_scalar_exch_coeff,
     1340!--          leaf_surface_conc from file which contains complete 3D data (separate vertical profiles
     1341!--          for each location).
    13381342             ELSE
    13391343                CALL pcm_read_plant_canopy_3d
    13401344             ENDIF
    13411345!
    1342 !--          Initialize LAD with data from file. If LAD is given in NetCDF file,
    1343 !--          use these values, else take LAD profiles from ASCII file.
    1344 !--          Please note, in NetCDF file LAD is only given up to the maximum
    1345 !--          canopy top, indicated by basal_area_density_f%nz. 
     1346!--          Initialize LAD with data from file. If LAD is given in NetCDF file, use these values,
     1347!--          else take LAD profiles from ASCII file.
     1348!--          Please note, in NetCDF file LAD is only given up to the maximum canopy top, indicated
     1349!--          by basal_area_density_f%nz.
    13461350             bad_s = 0.0_wp
    13471351             IF ( basal_area_density_f%from_file )  THEN
    13481352!
    1349 !--             Set also pch_index, used to be the upper bound of the vertical
    1350 !--             loops. Therefore, use the global top of the canopy layer.
     1353!--             Set also pch_index, used to be the upper bound of the vertical loops. Therefore, use
     1354!--             the global top of the canopy layer.
    13511355                pch_index = basal_area_density_f%nz - 1
    13521356
     
    13601364!--                   Check if resolved vegetation is mapped onto buildings.
    13611365!--                   Please see comment for leaf_area density
    1362                       IF ( ANY( bad_s(:,j,i) /= 0.0_wp )  .AND.                &
    1363                            ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND.  &
     1366                      IF ( ANY( bad_s(:,j,i) /= 0.0_wp )                 .AND.                     &
     1367                           ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) )  .AND.                     &
    13641368                           ANY( BTEST( wall_flags_total_0(:,j,i), 4 ) ) )  THEN
    13651369                         bad_s(:,j,i) = 0.0_wp
     
    13691373                ENDDO
    13701374#if defined( __parallel )
    1371                CALL MPI_ALLREDUCE( MPI_IN_PLACE, bad_on_top, 1, MPI_LOGICAL,  &
    1372                                    MPI_LOR, comm2d, ierr)
     1375               CALL MPI_ALLREDUCE( MPI_IN_PLACE, bad_on_top, 1, MPI_LOGICAL, MPI_LOR, comm2d, ierr)
    13731376#endif
    13741377                IF ( bad_on_top )  THEN
     
    13851388          CASE DEFAULT
    13861389!
    1387 !--          The DEFAULT case is reached either if the parameter
    1388 !--          canopy mode contains a wrong character string or if the
    1389 !--          user has coded a special case in the user interface.
    1390 !--          There, the subroutine user_init_plant_canopy checks
    1391 !--          which of these two conditions applies.
     1390!--          The DEFAULT case is reached either if the parameter canopy mode contains a wrong
     1391!--          character string or if the user has coded a special case in the user interface.
     1392!--          There, the subroutine user_init_plant_canopy checks which of these two conditions
     1393!--          applies.
    13921394             CALL user_init_plant_canopy
    1393  
     1395
    13941396       END SELECT
    13951397!
    1396 !--    Check that at least one grid point has an LAD /= 0, else this may
    1397 !--    cause errors in the radiation model.
     1398!--    Check that at least one grid point has an LAD /= 0, else this may cause errors in the
     1399!--    radiation model.
    13981400       lad_max = MAXVAL( lad_s )
    13991401#if defined( __parallel )
    1400        CALL MPI_ALLREDUCE( MPI_IN_PLACE, lad_max, 1, MPI_REAL, MPI_MAX,        &
    1401                            comm2d, ierr)
     1402       CALL MPI_ALLREDUCE( MPI_IN_PLACE, lad_max, 1, MPI_REAL, MPI_MAX, comm2d, ierr)
    14021403#endif
    14031404       IF ( lad_max <= 0.0_wp )  THEN
    1404           message_string = 'Plant-canopy model is switched-on but no ' //      &
     1405          message_string = 'Plant-canopy model is switched-on but no ' //                          &
    14051406                           'plant canopy is present in the model domain.'
    14061407          CALL message( 'pcm_init', 'PA0685', 1, 2, 0, 6, 0 )
    14071408       ENDIF
    1408    
    1409 !
    1410 !--    Initialize 2D index array indicating canopy top index. 
     1409
     1410!
     1411!--    Initialize 2D index array indicating canopy top index.
    14111412       ALLOCATE( pch_index_ji(nysg:nyng,nxlg:nxrg) )
    14121413       pch_index_ji = 0
    1413        
     1414
    14141415       DO  i = nxlg, nxrg
    14151416          DO  j = nysg, nyng
     
    14181419             ENDDO
    14191420!
    1420 !--          Check whether topography and local vegetation on top exceed
    1421 !--          height of the model domain.
     1421!--          Check whether topography and local vegetation on top exceed height of the model domain.
    14221422             IF ( topo_top_ind(j,i,0) + pch_index_ji(j,i) >= nzt + 1 )  THEN
    1423                 message_string =  'Local vegetation height on top of ' //      &
     1423                message_string =  'Local vegetation height on top of ' //                          &
    14241424                                  'topography exceeds height of model domain.'
    14251425                CALL message( 'pcm_init', 'PA0674', 2, 2, myid, 6, 0 )
     
    14341434!--    Exchange pch_index from all processors
    14351435#if defined( __parallel )
    1436        CALL MPI_ALLREDUCE( MPI_IN_PLACE, pch_index, 1, MPI_INTEGER,            &
    1437                            MPI_MAX, comm2d, ierr)
     1436       CALL MPI_ALLREDUCE( MPI_IN_PLACE, pch_index, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr)
    14381437#endif
    14391438!
     
    14411440       ALLOCATE( pcm_heating_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
    14421441       pcm_heating_rate = 0.0_wp
    1443        
     1442
    14441443       IF ( humidity )  THEN
    14451444          ALLOCATE( pcm_transpiration_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
     
    14491448       ENDIF
    14501449!
    1451 !--    Initialization of the canopy heat source distribution due to heating
    1452 !--    of the canopy layers by incoming solar radiation, in case that a non-zero
    1453 !--    value is set for the canopy top heat flux (cthf), which equals the
    1454 !--    available net radiation at canopy top.
    1455 !--    The heat source distribution is calculated by a decaying exponential
    1456 !--    function of the downward cumulative leaf area index (cum_lai_hf),
    1457 !--    assuming that the foliage inside the plant canopy is heated by solar
    1458 !--    radiation penetrating the canopy layers according to the distribution
    1459 !--    of net radiation as suggested by Brown & Covey (1966; Agric. Meteorol. 3,
    1460 !--    73–96). This approach has been applied e.g. by Shaw & Schumann (1992;
    1461 !--    Bound.-Layer Meteorol. 61, 47–64).
    1462 !--    When using the radiation_interactions, canopy heating (pcm_heating_rate)
    1463 !--    and plant canopy transpiration (pcm_transpiration_rate, pcm_latent_rate)
    1464 !--    are calculated in the RTM after the calculation of radiation.
     1450!--    Initialization of the canopy heat source distribution due to heating of the canopy layers by
     1451!--    incoming solar radiation, in case that a non-zero
     1452!--    value is set for the canopy top heat flux (cthf), which equals the available net radiation at
     1453!--    canopy top.
     1454!--    The heat source distribution is calculated by a decaying exponential function of the downward
     1455!--    cumulative leaf area index (cum_lai_hf), assuming that the foliage inside the plant canopy is
     1456!--    heated by solar radiation penetrating the canopy layers according to the distribution of net
     1457!--    radiation as suggested by Brown & Covey (1966; Agric. Meteorol. 3, 73–96). This approach has
     1458!--    been applied e.g. by Shaw & Schumann (1992; Bound.-Layer Meteorol. 61, 47–64).
     1459!--    When using the radiation_interactions, canopy heating (pcm_heating_rate) and plant canopy
     1460!--    transpiration (pcm_transpiration_rate, pcm_latent_rate) are calculated in the RTM after the
     1461!--    calculation of radiation.
    14651462       IF ( cthf /= 0.0_wp )  THEN
    14661463
    14671464          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    14681465!
    1469 !--       Piecewise calculation of the cumulative leaf area index by vertical
    1470 !--       integration of the leaf area density
     1466!--       Piecewise calculation of the cumulative leaf area index by vertical integration of the
     1467!--       leaf area density
    14711468          cum_lai_hf(:,:,:) = 0.0_wp
    14721469          DO  i = nxlg, nxrg
     
    14741471                DO  k = pch_index_ji(j,i)-1, 0, -1
    14751472                   IF ( k == pch_index_ji(j,i)-1 )  THEN
    1476                       cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
    1477                          ( 0.5_wp * lad_s(k+1,j,i) *                           &
    1478                            ( zw(k+1) - zu(k+1) ) )  +                          &
    1479                          ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
    1480                                                  lad_s(k,j,i) ) +              &
    1481                                       lad_s(k+1,j,i) ) *                       &
    1482                            ( zu(k+1) - zw(k) ) )
     1473                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                                    &
     1474                                          ( 0.5_wp * lad_s(k+1,j,i) *                              &
     1475                                            ( zw(k+1) - zu(k+1) ) )  +                             &
     1476                                          ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +                 &
     1477                                                                  lad_s(k,j,i) ) +                 &
     1478                                                       lad_s(k+1,j,i) ) *                          &
     1479                                            ( zu(k+1) - zw(k) ) )
    14831480                   ELSE
    1484                       cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
    1485                          ( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) +              &
    1486                                                  lad_s(k+1,j,i) ) +            &
    1487                                       lad_s(k+1,j,i) ) *                       &
    1488                            ( zw(k+1) - zu(k+1) ) )  +                          &
    1489                          ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
    1490                                                  lad_s(k,j,i) ) +              &
    1491                                       lad_s(k+1,j,i) ) *                       &
    1492                            ( zu(k+1) - zw(k) ) )
     1481                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                                    &
     1482                                          ( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) +                 &
     1483                                                                  lad_s(k+1,j,i) ) +               &
     1484                                                       lad_s(k+1,j,i) ) *                          &
     1485                                            ( zw(k+1) - zu(k+1) ) )  +                             &
     1486                                          ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +                 &
     1487                                                                  lad_s(k,j,i) ) +                 &
     1488                                                       lad_s(k+1,j,i) ) *                          &
     1489                                            ( zu(k+1) - zw(k) ) )
    14931490                   ENDIF
    14941491                ENDDO
     
    14971494
    14981495!
    1499 !--       In areas with canopy the surface value of the canopy heat
    1500 !--       flux distribution overrides the surface heat flux (shf)
     1496!--       In areas with canopy the surface value of the canopy heat flux distribution overrides the
     1497!--       surface heat flux (shf),
    15011498!--       Start with default surface type
    15021499          DO  m = 1, surf_def_h(0)%ns
    15031500             i = surf_def_h(0)%i(m)
    15041501             j = surf_def_h(0)%j(m)
    1505              IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
    1506                 surf_def_h(0)%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
     1502             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                                    &
     1503                surf_def_h(0)%shf(m) = cthf * EXP( -ext_coef * cum_lai_hf(0,j,i) )
    15071504          ENDDO
    15081505!
     
    15111508             i = surf_lsm_h(0)%i(m)
    15121509             j = surf_lsm_h(0)%j(m)
    1513              IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
    1514                 surf_lsm_h(0)%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
     1510             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                                    &
     1511                surf_lsm_h(0)%shf(m) = cthf * EXP( -ext_coef * cum_lai_hf(0,j,i) )
    15151512          ENDDO
    15161513!
     
    15191516             i = surf_usm_h(0)%i(m)
    15201517             j = surf_usm_h(0)%j(m)
    1521              IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
    1522                 surf_usm_h(0)%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
     1518             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                                    &
     1519                surf_usm_h(0)%shf(m) = cthf * EXP( -ext_coef * cum_lai_hf(0,j,i) )
    15231520          ENDDO
    15241521!
    15251522!
    1526 !--       Calculation of the heating rate (K/s) within the different layers of
    1527 !--       the plant canopy. Calculation is only necessary in areas covered with
    1528 !--       canopy.
    1529 !--       Within the different canopy layers the plant-canopy heating
    1530 !--       rate (pcm_heating_rate) is calculated as the vertical
    1531 !--       divergence of the canopy heat fluxes at the top and bottom
    1532 !--       of the respective layer
     1523!--       Calculation of the heating rate (K/s) within the different layers of the plant canopy.
     1524!--       Calculation is only necessary in areas covered with canopy.
     1525!--       Within the different canopy layers the plant-canopy heating rate (pcm_heating_rate) is
     1526!--       calculated as the vertical divergence of the canopy heat fluxes at the top and bottom of
     1527!--       the respective layer.
    15331528          DO  i = nxlg, nxrg
    15341529             DO  j = nysg, nyng
    15351530                DO  k = 1, pch_index_ji(j,i)
    15361531                   IF ( cum_lai_hf(0,j,i) /= 0.0_wp )  THEN
    1537                       pcm_heating_rate(k,j,i) = cthf *                          &
    1538                                 ( exp(-ext_coef*cum_lai_hf(k,j,i)) -           &
    1539                                   exp(-ext_coef*cum_lai_hf(k-1,j,i) ) ) / dzw(k)
     1532                      pcm_heating_rate(k,j,i) = cthf *                                             &
     1533                                                ( EXP( -ext_coef * cum_lai_hf(k,j,i) ) -           &
     1534                                                  EXP( -ext_coef * cum_lai_hf(k-1,j,i) ) ) / dzw(k)
    15401535                   ENDIF
    15411536                ENDDO
     
    15501545
    15511546
    1552 !------------------------------------------------------------------------------!
     1547!--------------------------------------------------------------------------------------------------!
    15531548! Description:
    15541549! ------------
    15551550!> Parin for &plant_canopy_parameters for plant canopy model
    1556 !------------------------------------------------------------------------------!
     1551!--------------------------------------------------------------------------------------------------!
    15571552    SUBROUTINE pcm_parin
    15581553
    1559        CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
    1560        
    1561        NAMELIST /plant_canopy_parameters/                                      &
    1562                                   alpha_lad, beta_lad, canopy_drag_coeff,      &
    1563                                   canopy_mode, cthf,                           &
    1564                                   lad_surface, lad_type_coef,                  &
    1565                                   lad_vertical_gradient,                       &
    1566                                   lad_vertical_gradient_level,                 &
    1567                                   lai_beta,                                    &
    1568                                   leaf_scalar_exch_coeff,                      &
    1569                                   leaf_surface_conc, pch_index,                &
    1570                                   plant_canopy_transpiration
    1571 
    1572        NAMELIST /canopy_par/      alpha_lad, beta_lad, canopy_drag_coeff,      &
    1573                                   canopy_mode, cthf,                           &
    1574                                   lad_surface, lad_type_coef,                  &
    1575                                   lad_vertical_gradient,                       &
    1576                                   lad_vertical_gradient_level,                 &
    1577                                   lai_beta,                                    &
    1578                                   leaf_scalar_exch_coeff,                      &
    1579                                   leaf_surface_conc, pch_index,                &
    1580                                   plant_canopy_transpiration
     1554       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
     1555
     1556       NAMELIST /plant_canopy_parameters/  alpha_lad,                                              &
     1557                                           beta_lad,                                               &
     1558                                           canopy_drag_coeff,                                      &
     1559                                           canopy_mode,                                            &
     1560                                           cthf,                                                   &
     1561                                           lad_surface,                                            &
     1562                                           lad_type_coef,                                          &
     1563                                           lad_vertical_gradient,                                  &
     1564                                           lad_vertical_gradient_level,                            &
     1565                                           lai_beta,                                               &
     1566                                           leaf_scalar_exch_coeff,                                 &
     1567                                           leaf_surface_conc,                                      &
     1568                                           pch_index,                                              &
     1569                                           plant_canopy_transpiration
     1570
     1571       NAMELIST /canopy_par/               alpha_lad,                                              &
     1572                                           beta_lad,                                               &
     1573                                           canopy_drag_coeff,                                      &
     1574                                           canopy_mode,                                            &
     1575                                           cthf,                                                   &
     1576                                           lad_surface,                                            &
     1577                                           lad_type_coef,                                          &
     1578                                           lad_vertical_gradient,                                  &
     1579                                           lad_vertical_gradient_level,                            &
     1580                                           lai_beta,                                               &
     1581                                           leaf_scalar_exch_coeff,                                 &
     1582                                           leaf_surface_conc,                                      &
     1583                                           pch_index,                                              &
     1584                                           plant_canopy_transpiration
    15811585
    15821586       line = ' '
     
    15841588!
    15851589!--    Try to find plant-canopy model package
    1586        REWIND ( 11 )
     1590       REWIND( 11 )
    15871591       line = ' '
    15881592       DO WHILE ( INDEX( line, '&plant_canopy_parameters' ) == 0 )
    1589           READ ( 11, '(A)', END=12 )  line
     1593          READ( 11, '(A)', END=12 )  line
    15901594       ENDDO
    1591        BACKSPACE ( 11 )
     1595       BACKSPACE( 11 )
    15921596
    15931597!
    15941598!--    Read user-defined namelist
    1595        READ ( 11, plant_canopy_parameters, ERR = 10 )
     1599       READ( 11, plant_canopy_parameters, ERR = 10 )
    15961600
    15971601!
     
    16061610!
    16071611!--    Try to find old namelist
    1608  12    REWIND ( 11 )
     1612 12    REWIND( 11 )
    16091613       line = ' '
    16101614       DO WHILE ( INDEX( line, '&canopy_par' ) == 0 )
    1611           READ ( 11, '(A)', END=14 )  line
     1615          READ( 11, '(A)', END=14 )  line
    16121616       ENDDO
    1613        BACKSPACE ( 11 )
     1617       BACKSPACE( 11 )
    16141618
    16151619!
    16161620!--    Read user-defined namelist
    1617        READ ( 11, canopy_par, ERR = 13, END = 14 )
    1618 
    1619        message_string = 'namelist canopy_par is deprecated and will be ' // &
    1620                      'removed in near future. Please use namelist ' //      &
    1621                      'plant_canopy_parameters instead'
     1621       READ( 11, canopy_par, ERR = 13, END = 14 )
     1622
     1623       message_string = 'namelist canopy_par is deprecated and will be ' //                        &
     1624                        'removed in near future. Please use namelist ' //                          &
     1625                        'plant_canopy_parameters instead'
    16221626       CALL message( 'pcm_parin', 'PA0487', 0, 1, 0, 6, 0 )
    16231627
     
    16371641
    16381642
    1639 !------------------------------------------------------------------------------!
     1643!--------------------------------------------------------------------------------------------------!
    16401644! Description:
    16411645! ------------
     
    16491653!> ...
    16501654!>
    1651 !> i.e. first line determines number of levels and further lines represent plant
    1652 !> canopy data, one line per column and variable. In each data line,
    1653 !> dtype represents variable to be set:
     1655!> i.e. first line determines number of levels and further lines represent plant canopy data, one
     1656!> line per column and variable. In each data line, dtype represents variable to be set:
    16541657!>
    16551658!> dtype=1: leaf area density (lad_s)
    16561659!> dtype=2....n: some additional plant canopy input data quantity
    16571660!>
    1658 !> Zeros are added automatically above num_levels until top of domain.  Any
    1659 !> non-specified (x,y) columns have zero values as default.
    1660 !------------------------------------------------------------------------------!
     1661!> Zeros are added automatically above num_levels until top of domain. Any non-specified (x,y)
     1662!> columns have zero values as default.
     1663!--------------------------------------------------------------------------------------------------!
    16611664    SUBROUTINE pcm_read_plant_canopy_3d
    16621665
    1663        USE exchange_horiz_mod,                                                    &
     1666       USE exchange_horiz_mod,                                                                     &
    16641667           ONLY:  exchange_horiz
    16651668
    1666        INTEGER(iwp)                        ::  dtype     !< type of input data (1=lad)
    1667        INTEGER(iwp)                        ::  pctype    !< type of plant canopy (deciduous,non-deciduous,...)
    1668        INTEGER(iwp)                        ::  i, j      !< running index
    1669        INTEGER(iwp)                        ::  nzp       !< number of vertical layers of plant canopy
    1670        INTEGER(iwp)                        ::  nzpltop   !<
    1671        INTEGER(iwp)                        ::  nzpl      !<
    1672        INTEGER(iwp)                        ::  kk        !<
    1673        
     1669       INTEGER(iwp) ::  dtype     !< type of input data (1=lad)
     1670       INTEGER(iwp) ::  i         !< running index
     1671       INTEGER(iwp) ::  j         !< running index
     1672       INTEGER(iwp) ::  kk        !<
     1673       INTEGER(iwp) ::  nzp       !< number of vertical layers of plant canopy
     1674       INTEGER(iwp) ::  nzpltop   !<
     1675       INTEGER(iwp) ::  nzpl      !<
     1676       INTEGER(iwp) ::  pctype    !< type of plant canopy (deciduous,non-deciduous,...)
     1677
    16741678       REAL(wp), DIMENSION(:), ALLOCATABLE ::  col   !< vertical column of input data
    16751679
     
    16771681!--    Initialize lad_s array
    16781682       lad_s = 0.0_wp
    1679        
     1683
    16801684!
    16811685!--    Open and read plant canopy input data
    1682        OPEN(152, FILE='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ),         &
    1683                  ACCESS='SEQUENTIAL', ACTION='READ', STATUS='OLD',             &
    1684                  FORM='FORMATTED', ERR=515)
    1685        READ(152, *, ERR=516, END=517)  nzp   !< read first line = number of vertical layers
     1686       OPEN( 152, FILE='PLANT_CANOPY_DATA_3D' // TRIM( coupling_char ), ACCESS='SEQUENTIAL',       &
     1687                  ACTION='READ', STATUS='OLD', FORM='FORMATTED', ERR=515 )
     1688       READ( 152, *, ERR=516, END=517 )  nzp   !< read first line = number of vertical layers
    16861689       nzpltop = MIN(nzt+1, nzb+nzp-1)
    16871690       nzpl = nzpltop - nzb + 1    !< no. of layers to assign
     
    16891692
    16901693       DO
    1691           READ(152, *, ERR=516, END=517) dtype, i, j, pctype, col(:)
     1694          READ( 152, *, ERR=516, END=517 ) dtype, i, j, pctype, col(:)
    16921695          IF ( i < nxlg  .OR.  i > nxrg  .OR.  j < nysg  .OR.  j > nyng )  CYCLE
    1693          
    1694           SELECT CASE (dtype)
     1696
     1697          SELECT CASE ( dtype )
    16951698             CASE( 1 )   !< leaf area density
    16961699!
    1697 !--             This is just the pure canopy layer assumed to be grounded to
    1698 !--             a flat domain surface. At locations where plant canopy sits
    1699 !--             on top of any kind of topography, the vertical plant column
    1700 !--             must be "lifted", which is done in SUBROUTINE pcm_tendency.           
     1700!--             This is just the pure canopy layer assumed to be grounded to a flat domain surface.
     1701!--             At locations where plant canopy sits on top of any kind of topography, the vertical
     1702!--             plant column must be "lifted", which is done in SUBROUTINE pcm_tendency.
    17011703                IF ( pctype < 0  .OR.  pctype > 10 )  THEN   !< incorrect plant canopy type
    1702                    WRITE( message_string, * ) 'Incorrect type of plant canopy. '   //  &
    1703                                               'Allowed values 0 <= pctype <= 10, ' //  &
     1704                   WRITE( message_string, * ) 'Incorrect type of plant canopy. '   //              &
     1705                                              'Allowed values 0 <= pctype <= 10, ' //              &
    17041706                                              'but pctype is ', pctype
    17051707                   CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 )
     
    17081710                lad_s(nzb:nzpltop-kk, j, i) = col(kk:nzpl-1)*lad_type_coef(pctype)
    17091711             CASE DEFAULT
    1710                 WRITE(message_string, '(a,i2,a)')   &
     1712                WRITE( message_string, '(a,i2,a)' )                                                &
    17111713                     'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
    17121714                CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
     
    17201722       CALL message( 'pcm_read_plant_canopy_3d', 'PA0532', 1, 2, 0, 6, 0 )
    17211723
    1722 517    CLOSE(152)
     1724517    CLOSE( 152 )
    17231725       DEALLOCATE( col )
    1724        
     1726
    17251727       CALL exchange_horiz( lad_s, nbgp )
    1726        
     1728
    17271729    END SUBROUTINE pcm_read_plant_canopy_3d
    17281730
    1729 !------------------------------------------------------------------------------!
     1731!--------------------------------------------------------------------------------------------------!
    17301732! Description:
    17311733! ------------
    17321734!> Read module-specific global restart data (Fortran binary format).
    1733 !------------------------------------------------------------------------------!
     1735!--------------------------------------------------------------------------------------------------!
    17341736    SUBROUTINE pcm_rrd_global_ftn( found )
    17351737
    1736        LOGICAL, INTENT(OUT)  ::  found
     1738       LOGICAL, INTENT(OUT) ::  found
     1739
    17371740
    17381741       found = .TRUE.
     
    17411744
    17421745          CASE ( 'pch_index' )
    1743              READ ( 13 )  pch_index
     1746             READ( 13 )  pch_index
    17441747
    17451748          CASE DEFAULT
     
    17511754    END SUBROUTINE pcm_rrd_global_ftn
    17521755
    1753 !------------------------------------------------------------------------------!
     1756!--------------------------------------------------------------------------------------------------!
    17541757! Description:
    17551758! ------------
    17561759!> Read module-specific global restart data (MPI-IO).
    1757 !------------------------------------------------------------------------------!
     1760!--------------------------------------------------------------------------------------------------!
    17581761    SUBROUTINE pcm_rrd_global_mpi
    17591762
     
    17631766
    17641767
    1765 !------------------------------------------------------------------------------!
     1768!--------------------------------------------------------------------------------------------------!
    17661769! Description:
    17671770! ------------
    17681771!> Read module-specific local restart data arrays (Fortran binary format).
    1769 !------------------------------------------------------------------------------!
    1770     SUBROUTINE pcm_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,          &
    1771                                   nxr_on_file, nynf, nync, nyn_on_file, nysf,      &
    1772                                   nysc, nys_on_file, found )
    1773 
    1774        INTEGER(iwp) ::  k               !<
    1775        INTEGER(iwp) ::  nxlc            !<
    1776        INTEGER(iwp) ::  nxlf            !<
    1777        INTEGER(iwp) ::  nxl_on_file     !<
    1778        INTEGER(iwp) ::  nxrc            !<
    1779        INTEGER(iwp) ::  nxrf            !<
    1780        INTEGER(iwp) ::  nxr_on_file     !<
    1781        INTEGER(iwp) ::  nync            !<
    1782        INTEGER(iwp) ::  nynf            !<
    1783        INTEGER(iwp) ::  nyn_on_file     !<
    1784        INTEGER(iwp) ::  nysc            !<
    1785        INTEGER(iwp) ::  nysf            !<
    1786        INTEGER(iwp) ::  nys_on_file     !<
    1787 
    1788        LOGICAL, INTENT(OUT)  :: found
    1789 
    1790        REAL(wp), DIMENSION(0:pch_index,                                        &
    1791                            nys_on_file-nbgp:nyn_on_file+nbgp,                  &
    1792                            nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d2  !< temporary 3D array for entire vertical
    1793                                                                           !< extension of canopy layer
     1772!--------------------------------------------------------------------------------------------------!
     1773    SUBROUTINE pcm_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync, &
     1774                                  nyn_on_file, nysf, nysc, nys_on_file, found )
     1775
     1776       INTEGER(iwp) ::  k               !<
     1777       INTEGER(iwp) ::  nxl_on_file     !<
     1778       INTEGER(iwp) ::  nxlc            !<
     1779       INTEGER(iwp) ::  nxlf            !<
     1780       INTEGER(iwp) ::  nxr_on_file     !<
     1781       INTEGER(iwp) ::  nxrc            !<
     1782       INTEGER(iwp) ::  nxrf            !<
     1783       INTEGER(iwp) ::  nyn_on_file     !<
     1784       INTEGER(iwp) ::  nync            !<
     1785       INTEGER(iwp) ::  nynf            !<
     1786       INTEGER(iwp) ::  nys_on_file     !<
     1787       INTEGER(iwp) ::  nysc            !<
     1788       INTEGER(iwp) ::  nysf            !<
     1789
     1790       LOGICAL, INTENT(OUT) ::  found
     1791
     1792       REAL(wp), DIMENSION( 0:pch_index,                                                           &
     1793                            nys_on_file-nbgp:nyn_on_file+nbgp,                                     &
     1794                            nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d2  !< temporary 3D array for entire vertical
     1795                                                                           !< extension of canopy layer
    17941796       found = .TRUE.
     1797
    17951798
    17961799       SELECT CASE ( restart_string(1:length) )
     
    18001803                ALLOCATE( pcm_heatrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
    18011804                pcm_heatrate_av = 0.0_wp
    1802              ENDIF 
    1803              IF ( k == 1 )  READ ( 13 )  tmp_3d2
    1804              pcm_heatrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     1805             ENDIF
     1806             IF ( k == 1 )  READ( 13 )  tmp_3d2
     1807             pcm_heatrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                &
    18051808                        tmp_3d2(0:pch_index,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    18061809
     
    18091812                ALLOCATE( pcm_latentrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
    18101813                pcm_latentrate_av = 0.0_wp
    1811              ENDIF 
    1812              IF ( k == 1 )  READ ( 13 )  tmp_3d2
    1813              pcm_latentrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     1814             ENDIF
     1815             IF ( k == 1 )  READ( 13 )  tmp_3d2
     1816             pcm_latentrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =              &
    18141817                        tmp_3d2(0:pch_index,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    18151818
     
    18181821                ALLOCATE( pcm_transpirationrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
    18191822                pcm_transpirationrate_av = 0.0_wp
    1820              ENDIF 
    1821              IF ( k == 1 )  READ ( 13 )  tmp_3d2
    1822              pcm_transpirationrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     1823             ENDIF
     1824             IF ( k == 1 )  READ( 13 )  tmp_3d2
     1825             pcm_transpirationrate_av(0:pch_index,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =       &
    18231826                        tmp_3d2(0:pch_index,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    18241827
     
    18321835
    18331836
    1834 !------------------------------------------------------------------------------!
     1837!--------------------------------------------------------------------------------------------------!
    18351838! Description:
    18361839! ------------
    18371840!> Read module-specific local restart data arrays (MPI-IO).
    1838 !------------------------------------------------------------------------------!
     1841!--------------------------------------------------------------------------------------------------!
    18391842    SUBROUTINE pcm_rrd_local_mpi
    18401843
     
    18851888
    18861889
    1887 !------------------------------------------------------------------------------!
     1890!--------------------------------------------------------------------------------------------------!
    18881891! Description:
    18891892! ------------
    1890 !> Calculation of the tendency terms, accounting for the effect of the plant
    1891 !> canopy on momentum and scalar quantities.
     1893!> Calculation of the tendency terms, accounting for the effect of the plant canopy on momentum and
     1894!> scalar quantities.
    18921895!>
    1893 !> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
    1894 !> (defined on scalar grid), as initialized in subroutine pcm_init.
    1895 !> The lad on the w-grid is vertically interpolated from the surrounding
    1896 !> lad_s. The upper boundary of the canopy is defined on the w-grid at
    1897 !> k = pch_index. Here, the lad is zero.
     1896!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0 (defined on scalar grid), as
     1897!> initialized in subroutine pcm_init.
     1898!> The lad on the w-grid is vertically interpolated from the surrounding lad_s. The upper boundary
     1899!> of the canopy is defined on the w-grid at k = pch_index. Here, the lad is zero.
    18981900!>
    1899 !> The canopy drag must be limited (previously accounted for by calculation of
    1900 !> a limiting canopy timestep for the determination of the maximum LES timestep
    1901 !> in subroutine timestep), since it is physically impossible that the canopy
    1902 !> drag alone can locally change the sign of a velocity component. This
    1903 !> limitation is realized by calculating preliminary tendencies and velocities.
    1904 !> It is subsequently checked if the preliminary new velocity has a different
    1905 !> sign than the current velocity. If so, the tendency is limited in a way that
    1906 !> the velocity can at maximum be reduced to zero by the canopy drag.
     1901!> The canopy drag must be limited (previously accounted for by calculation of a limiting canopy
     1902!> timestep for the determination of the maximum LES timestep in subroutine timestep), since it is
     1903!> physically impossible that the canopy drag alone can locally change the sign of a velocity
     1904!> component. This limitation is realized by calculating preliminary tendencies and velocities. It
     1905!> is subsequently checked if the preliminary new velocity has a different sign than the current
     1906!> velocity. If so, the tendency is limited in a way that the velocity can at maximum be reduced to
     1907!> zero by the canopy drag.
    19071908!>
    19081909!>
    19091910!> Call for all grid points
    1910 !------------------------------------------------------------------------------!
     1911!--------------------------------------------------------------------------------------------------!
    19111912    SUBROUTINE pcm_tendency( component )
    19121913
     
    19261927       REAL(wp) ::  lad_local !< local lad value
    19271928       REAL(wp) ::  pre_tend  !< preliminary tendency
    1928        REAL(wp) ::  pre_u     !< preliminary u-value 
    1929        REAL(wp) ::  pre_v     !< preliminary v-value 
    1930        REAL(wp) ::  pre_w     !< preliminary w-value 
     1929       REAL(wp) ::  pre_u     !< preliminary u-value
     1930       REAL(wp) ::  pre_v     !< preliminary v-value
     1931       REAL(wp) ::  pre_w     !< preliminary w-value
    19311932
    19321933
    19331934       ddt_3d = 1.0_wp / dt_3d
    1934  
     1935
    19351936!
    19361937!--    Compute drag for the three velocity components and the SGS-TKE:
     
    19431944                DO  j = nys, nyn
    19441945!
    1945 !--                Set control flags indicating east- and westward-orientated
    1946 !--                building edges. Note, building_egde_w is set from the perspective
    1947 !--                of the potential rooftop grid point, while building_edge_e is
    1948 !--                is set from the perspective of the non-building grid point.
    1949                    building_edge_w = ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )&
    1950                         .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )
    1951                    building_edge_e = ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )&
    1952                         .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
     1946!--                Set control flags indicating east- and westward-orientated building edges. Note,
     1947!--                building_egde_w is set from the perspective of the potential rooftop grid point,
     1948!--                while building_edge_e is set from the perspective of the non-building grid point.
     1949                   building_edge_w =              ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )   &
     1950                                     .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )
     1951                   building_edge_e =              ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )   &
     1952                                     .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
    19531953!
    19541954!--                Determine topography-top index on u-grid
     
    19571957                      kk = k - topo_top_ind(j,i,1)   !- lad arrays are defined flat
    19581958!
    1959 !--                   In order to create sharp boundaries of the plant canopy,
    1960 !--                   the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
    1961 !--                   rather than being interpolated from the surrounding lad_s,
    1962 !--                   because this would yield smaller lad at the canopy boundaries
    1963 !--                   than inside of the canopy.
    1964 !--                   For the same reason, the lad at the rightmost(i+1)canopy
    1965 !--                   boundary on the u-grid equals lad_s(k,j,i), which is considered
    1966 !--                   in the next if-statement. Note, at left-sided building edges
    1967 !--                   this is not applied, here the LAD is equals the LAD at grid
    1968 !--                   point (k,j,i), in order to avoid that LAD is mistakenly mapped
    1969 !--                   on top of a roof where (usually) is no LAD is defined. The same is
    1970 !--                   also valid for bad_s.
     1959!--                   In order to create sharp boundaries of the plant canopy, the lad on the u-grid
     1960!--                   at index (k,j,i) is equal to lad_s(k,j,i), rather than being interpolated from
     1961!--                   the surrounding lad_s, because this would yield smaller lad at the canopy
     1962!--                   boundaries than inside of the canopy.
     1963!--                   For the same reason, the lad at the rightmost(i+1)canopy boundary on the
     1964!--                   u-grid equals lad_s(k,j,i), which is considered in the next if-statement.
     1965!--                   Note, at left-sided building edges this is not applied, here the LAD equals
     1966!--                   the LAD at grid point (k,j,i), in order to avoid that LAD is mistakenly mapped
     1967!--                   on top of a roof where (usually) no LAD is defined. The same is also valid for
     1968!--                   bad_s.
    19711969                      lad_local = lad_s(kk,j,i)
    19721970                      IF ( lad_local == 0.0_wp  .AND.  lad_s(kk,j,i-1) > 0.0_wp                    &
     
    19771975                           .AND.  .NOT. building_edge_w )  bad_local = bad_s(kk,j,i-1)
    19781976!
    1979 !--                   In order to avoid that LAD is mistakenly considered at right-
    1980 !--                   sided building edges (here the topography-top index for the
    1981 !--                   u-component at index j,i is still on the building while the
    1982 !--                   topography top for the scalar isn't), LAD is taken from grid
    1983 !--                   point (j,i-1). The same is also valid for bad_s.
     1977!--                   In order to avoid that LAD is mistakenly considered at right-sided building
     1978!--                   edges (here the topography-top index for the u-component at index j,i is still
     1979!--                   on the building while the topography top for the scalar isn't), LAD is taken
     1980!--                   from grid point (j,i-1). The same is also valid for bad_s.
    19841981                      IF ( lad_local > 0.0_wp  .AND.  lad_s(kk,j,i-1) == 0.0_wp                    &
    19851982                           .AND.  building_edge_e )  lad_local = lad_s(kk,j,i-1)
     
    19911988!
    19921989!--                   Calculate preliminary value (pre_tend) of the tendency
    1993                       pre_tend = - canopy_drag_coeff *                         &
    1994                                    ( lad_local + bad_local ) *                 &
    1995                                    SQRT( u(k,j,i)**2 +                         &
    1996                                          ( 0.25_wp * ( v(k,j,i-1) +            &
    1997                                                        v(k,j,i)   +            &
    1998                                                        v(k,j+1,i) +            &
    1999                                                        v(k,j+1,i-1) )          &
    2000                                          )**2 +                                &
    2001                                          ( 0.25_wp * ( w(k-1,j,i-1) +          &
    2002                                                        w(k-1,j,i)   +          &
    2003                                                        w(k,j,i-1)   +          &
    2004                                                        w(k,j,i) )              &
    2005                                          )**2                                  &
    2006                                        ) *                                     &
     1990                      pre_tend = - canopy_drag_coeff *                                             &
     1991                                   ( lad_local + bad_local ) *                                     &
     1992                                   SQRT( u(k,j,i)**2 +                                             &
     1993                                         ( 0.25_wp * ( v(k,j,i-1) +                                &
     1994                                                       v(k,j,i)   +                                &
     1995                                                       v(k,j+1,i) +                                &
     1996                                                       v(k,j+1,i-1) )                              &
     1997                                         )**2 +                                                    &
     1998                                         ( 0.25_wp * ( w(k-1,j,i-1) +                              &
     1999                                                       w(k-1,j,i)   +                              &
     2000                                                       w(k,j,i-1)   +                              &
     2001                                                       w(k,j,i) )                                  &
     2002                                         )**2                                                      &
     2003                                       ) *                                                         &
    20072004                                   u(k,j,i)
    20082005
     
    20132010!--                   Compare sign of old velocity and new preliminary velocity,
    20142011!--                   and in case the signs are different, limit the tendency
    2015                       IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
     2012                      IF ( SIGN( pre_u,u(k,j,i) ) /= pre_u )  THEN
    20162013                         pre_tend = - u(k,j,i) * ddt_3d
    20172014                      ENDIF
     
    20302027                DO  j = nysv, nyn
    20312028!
    2032 !--                Set control flags indicating north- and southward-orientated
    2033 !--                building edges. Note, building_egde_s is set from the perspective
    2034 !--                of the potential rooftop grid point, while building_edge_n is
    2035 !--                is set from the perspective of the non-building grid point.
    2036                    building_edge_s = ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )&
    2037                         .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )
    2038                    building_edge_n = ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )&
    2039                         .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
     2029!--                Set control flags indicating north- and southward-orientated building edges.
     2030!--                Note, building_egde_s is set from the perspective of the potential rooftop grid
     2031!--                point, while building_edge_n is set from the perspective of the non-building grid
     2032!--                point.
     2033                   building_edge_s =              ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )   &
     2034                                     .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )
     2035                   building_edge_n =              ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )   &
     2036                                     .AND.  .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
    20402037!
    20412038!--                Determine topography-top index on v-grid
     
    20442041                      kk = k - topo_top_ind(j,i,2)   !- lad arrays are defined flat
    20452042!
    2046 !--                   In order to create sharp boundaries of the plant canopy,
    2047 !--                   the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
    2048 !--                   rather than being interpolated from the surrounding lad_s,
    2049 !--                   because this would yield smaller lad at the canopy boundaries
    2050 !--                   than inside of the canopy.
    2051 !--                   For the same reason, the lad at the northmost(j+1)canopy
    2052 !--                   boundary on the v-grid equals lad_s(k,j,i), which is considered
    2053 !--                   in the next if-statement. Note, at left-sided building edges
    2054 !--                   this is not applied, here the LAD is equals the LAD at grid
    2055 !--                   point (k,j,i), in order to avoid that LAD is mistakenly mapped
    2056 !--                   on top of a roof where (usually) is no LAD is defined.
     2043!--                   In order to create sharp boundaries of the plant canopy, the lad on the v-grid
     2044!--                   at index (k,j,i) is equal to lad_s(k,j,i), rather than being interpolated from
     2045!--                   the surrounding lad_s, because this would yield smaller lad at the canopy
     2046!--                   boundaries  than inside of the canopy.
     2047!--                   For the same reason, the lad at the northmost (j+1) canopy boundary on the
     2048!--                   v-grid equals lad_s(k,j,i), which is considered in the next if-statement.
     2049!--                   Note, at left-sided building edges this is not applied, here the LAD equals
     2050!--                   the LAD at grid point (k,j,i), in order to avoid that LAD is mistakenly mapped
     2051!--                   on top of a roof where (usually) no LAD is defined.
    20572052!--                   The same is also valid for bad_s.
    20582053                      lad_local = lad_s(kk,j,i)
    20592054                      IF ( lad_local == 0.0_wp  .AND.  lad_s(kk,j-1,i) > 0.0_wp                    &
    2060                        .AND.  .NOT. building_edge_s )  lad_local = lad_s(kk,j-1,i)
     2055                           .AND.  .NOT. building_edge_s )  lad_local = lad_s(kk,j-1,i)
    20612056
    20622057                      bad_local = bad_s(kk,j,i)
    20632058                      IF ( bad_local == 0.0_wp  .AND.  bad_s(kk,j-1,i) > 0.0_wp                    &
    2064                        .AND.  .NOT. building_edge_s )  bad_local = bad_s(kk,j-1,i)
    2065 !
    2066 !--                   In order to avoid that LAD is mistakenly considered at right-
    2067 !--                   sided building edges (here the topography-top index for the
    2068 !--                   u-component at index j,i is still on the building while the
    2069 !--                   topography top for the scalar isn't), LAD is taken from grid
    2070 !--                   point (j,i-1). The same is also valid for bad_s.
     2059                           .AND.  .NOT. building_edge_s )  bad_local = bad_s(kk,j-1,i)
     2060!
     2061!--                   In order to avoid that LAD is mistakenly considered at right-sided building
     2062!--                   edges (here the topography-top index for the u-component at index j,i is still
     2063!--                   on the building while the topography top for the scalar isn't), LAD is taken
     2064!--                   from grid point (j,i-1). The same is also valid for bad_s.
    20712065                      IF ( lad_local > 0.0_wp  .AND.  lad_s(kk,j-1,i) == 0.0_wp                    &
    2072                        .AND.  building_edge_n )  lad_local = lad_s(kk,j-1,i)
     2066                           .AND.  building_edge_n )  lad_local = lad_s(kk,j-1,i)
    20732067
    20742068                      IF ( bad_local > 0.0_wp  .AND.  bad_s(kk,j-1,i) == 0.0_wp                    &
    2075                        .AND.  building_edge_n )  bad_local = bad_s(kk,j-1,i)
     2069                           .AND.  building_edge_n )  bad_local = bad_s(kk,j-1,i)
    20762070
    20772071                      pre_tend = 0.0_wp
     
    20792073!
    20802074!--                   Calculate preliminary value (pre_tend) of the tendency
    2081                       pre_tend = - canopy_drag_coeff *                         &
    2082                                    ( lad_local + bad_local ) *                 &
    2083                                    SQRT( ( 0.25_wp * ( u(k,j-1,i)   +          &
    2084                                                        u(k,j-1,i+1) +          &
    2085                                                        u(k,j,i)     +          &
    2086                                                        u(k,j,i+1) )            &
    2087                                          )**2 +                                &
    2088                                          v(k,j,i)**2 +                         &
    2089                                          ( 0.25_wp * ( w(k-1,j-1,i) +          &
    2090                                                        w(k-1,j,i)   +          &
    2091                                                        w(k,j-1,i)   +          &
    2092                                                        w(k,j,i) )              &
    2093                                          )**2                                  &
    2094                                        ) *                                     &
     2075                      pre_tend = - canopy_drag_coeff *                                             &
     2076                                   ( lad_local + bad_local ) *                                     &
     2077                                   SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                              &
     2078                                                       u(k,j-1,i+1) +                              &
     2079                                                       u(k,j,i)     +                              &
     2080                                                       u(k,j,i+1) )                                &
     2081                                         )**2 +                                                    &
     2082                                         v(k,j,i)**2 +                                             &
     2083                                         ( 0.25_wp * ( w(k-1,j-1,i) +                              &
     2084                                                       w(k-1,j,i)   +                              &
     2085                                                       w(k,j-1,i)   +                              &
     2086                                                       w(k,j,i) )                                  &
     2087                                         )**2                                                      &
     2088                                       ) *                                                         &
    20952089                                   v(k,j,i)
    20962090
     
    20992093                      pre_v = v(k,j,i) + dt_3d * pre_tend
    21002094!
    2101 !--                   Compare sign of old velocity and new preliminary velocity,
    2102 !--                   and in case the signs are different, limit the tendency
    2103                       IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
     2095!--                   Compare sign of old velocity and new preliminary velocity, and in case the
     2096!--                   signs are different, limit the tendency.
     2097                      IF ( SIGN( pre_v,v(k,j,i) ) /= pre_v )  THEN
    21042098                         pre_tend = - v(k,j,i) * ddt_3d
    21052099                      ELSE
     
    21292123!
    21302124!--                   Calculate preliminary value (pre_tend) of the tendency
    2131                       pre_tend = - canopy_drag_coeff *                         &
    2132                                    ( 0.5_wp *                                  &
    2133                                       ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) +    &
    2134                                      0.5_wp *                                  &
    2135                                       ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) *  &
    2136                                    SQRT( ( 0.25_wp * ( u(k,j,i)   +            &
    2137                                                        u(k,j,i+1) +            &
    2138                                                        u(k+1,j,i) +            &
    2139                                                        u(k+1,j,i+1) )          &
    2140                                          )**2 +                                &
    2141                                          ( 0.25_wp * ( v(k,j,i)   +            &
    2142                                                        v(k,j+1,i) +            &
    2143                                                        v(k+1,j,i) +            &
    2144                                                        v(k+1,j+1,i) )          &
    2145                                          )**2 +                                &
    2146                                          w(k,j,i)**2                           &
    2147                                        ) *                                     &
     2125                      pre_tend = - canopy_drag_coeff *                                             &
     2126                                   ( 0.5_wp * ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) +                &
     2127                                     0.5_wp * ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) *              &
     2128                                   SQRT( ( 0.25_wp * ( u(k,j,i)   +                                &
     2129                                                       u(k,j,i+1) +                                &
     2130                                                       u(k+1,j,i) +                                &
     2131                                                       u(k+1,j,i+1) )                              &
     2132                                         )**2 +                                                    &
     2133                                         ( 0.25_wp * ( v(k,j,i)   +                                &
     2134                                                       v(k,j+1,i) +                                &
     2135                                                       v(k+1,j,i) +                                &
     2136                                                       v(k+1,j+1,i) )                              &
     2137                                         )**2 +                                                    &
     2138                                         w(k,j,i)**2                                               &
     2139                                       ) *                                                         &
    21482140                                   w(k,j,i)
    21492141!
     
    21512143                      pre_w = w(k,j,i) + dt_3d * pre_tend
    21522144!
    2153 !--                   Compare sign of old velocity and new preliminary velocity,
    2154 !--                   and in case the signs are different, limit the tendency
    2155                       IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
     2145!--                   Compare sign of old velocity and new preliminary velocity, and in case the
     2146!--                   signs are different, limit the tendency
     2147                      IF ( SIGN( pre_w,w(k,j,i) ) /= pre_w )  THEN
    21562148                         pre_tend = - w(k,j,i) * ddt_3d
    21572149                      ELSE
     
    21692161!--       potential temperature
    21702162          CASE ( 4 )
    2171              IF ( humidity ) THEN
     2163             IF ( humidity )  THEN
    21722164                DO  i = nxl, nxr
    21732165                   DO  j = nys, nyn
     
    22022194                      kk = k - topo_top_ind(j,i,0)   !- lad arrays are defined flat
    22032195
    2204                       IF ( .NOT. plant_canopy_transpiration ) THEN
     2196                      IF ( .NOT. plant_canopy_transpiration )  THEN
    22052197                         ! pcm_transpiration_rate is calculated in radiation model
    22062198                         ! in case of plant_canopy_transpiration = .T.
    22072199                         ! to include also the dependecy to the radiation
    22082200                         ! in the plant canopy box
    2209                          pcm_transpiration_rate(kk,j,i) =  - leaf_scalar_exch_coeff &
    2210                                           * lad_s(kk,j,i) *                         &
    2211                                           SQRT( ( 0.5_wp * ( u(k,j,i) +             &
    2212                                                              u(k,j,i+1) )           &
    2213                                                 )**2 +                              &
    2214                                                 ( 0.5_wp * ( v(k,j,i) +             &
    2215                                                              v(k,j+1,i) )           &
    2216                                                 )**2 +                              &
    2217                                                 ( 0.5_wp * ( w(k-1,j,i) +           &
    2218                                                              w(k,j,i) )             &
    2219                                                 )**2                                &
    2220                                               ) *                                   &
    2221                                           ( q(k,j,i) - leaf_surface_conc )
     2201                         pcm_transpiration_rate(kk,j,i) = - leaf_scalar_exch_coeff                &
     2202                                                            * lad_s(kk,j,i) *                      &
     2203                                                            SQRT( ( 0.5_wp * ( u(k,j,i) +          &
     2204                                                                               u(k,j,i+1) )        &
     2205                                                                  )**2 +                           &
     2206                                                                  ( 0.5_wp * ( v(k,j,i) +          &
     2207                                                                               v(k,j+1,i) )        &
     2208                                                                  )**2 +                           &
     2209                                                                  ( 0.5_wp * ( w(k-1,j,i) +        &
     2210                                                                               w(k,j,i) )          &
     2211                                                                  )**2                             &
     2212                                                                ) *                                &
     2213                                                            ( q(k,j,i) - leaf_surface_conc )
    22222214                      ENDIF
    22232215
     
    22372229
    22382230                      kk = k - topo_top_ind(j,i,0)   !- lad arrays are defined flat
    2239                       tend(k,j,i) = tend(k,j,i) -                              &
    2240                                        2.0_wp * canopy_drag_coeff *            &
    2241                                        ( lad_s(kk,j,i) + bad_s(kk,j,i) ) *     &
    2242                                        SQRT( ( 0.5_wp * ( u(k,j,i) +           &
    2243                                                           u(k,j,i+1) )         &
    2244                                              )**2 +                            &
    2245                                              ( 0.5_wp * ( v(k,j,i) +           &
    2246                                                           v(k,j+1,i) )         &
    2247                                              )**2 +                            &
    2248                                              ( 0.5_wp * ( w(k,j,i) +           &
    2249                                                           w(k+1,j,i) )         &
    2250                                              )**2                              &
    2251                                            ) *                                 &
    2252                                        e(k,j,i)
     2231                      tend(k,j,i) = tend(k,j,i) -  2.0_wp * canopy_drag_coeff *                    &
     2232                                                   ( lad_s(kk,j,i) + bad_s(kk,j,i) ) *             &
     2233                                                   SQRT( ( 0.5_wp * ( u(k,j,i) +                   &
     2234                                                                      u(k,j,i+1) )                 &
     2235                                                         )**2 +                                    &
     2236                                                         ( 0.5_wp * ( v(k,j,i) +                   &
     2237                                                                      v(k,j+1,i) )                 &
     2238                                                         )**2 +                                    &
     2239                                                         ( 0.5_wp * ( w(k,j,i) +                   &
     2240                                                                      w(k+1,j,i) )                 &
     2241                                                         )**2                                      &
     2242                                                       ) *                                         &
     2243                                                   e(k,j,i)
    22532244                   ENDDO
    22542245                ENDDO
    2255              ENDDO 
     2246             ENDDO
    22562247!
    22572248!--       scalar concentration
     
    22642255
    22652256                      kk = k - topo_top_ind(j,i,0)   !- lad arrays are defined flat
    2266                       tend(k,j,i) = tend(k,j,i) -                              &
    2267                                        leaf_scalar_exch_coeff *                &
    2268                                        lad_s(kk,j,i) *                         &
    2269                                        SQRT( ( 0.5_wp * ( u(k,j,i) +           &
    2270                                                           u(k,j,i+1) )         &
    2271                                              )**2 +                            &
    2272                                              ( 0.5_wp * ( v(k,j,i) +           &
    2273                                                           v(k,j+1,i) )         &
    2274                                              )**2 +                            &
    2275                                              ( 0.5_wp * ( w(k-1,j,i) +         &
    2276                                                           w(k,j,i) )           &
    2277                                              )**2                              &
    2278                                            ) *                                 &
    2279                                        ( s(k,j,i) - leaf_surface_conc )
     2257                      tend(k,j,i) = tend(k,j,i) -  leaf_scalar_exch_coeff *                        &
     2258                                                   lad_s(kk,j,i) *                                 &
     2259                                                   SQRT( ( 0.5_wp * ( u(k,j,i) +                   &
     2260                                                                      u(k,j,i+1) )                 &
     2261                                                         )**2 +                                    &
     2262                                                         ( 0.5_wp * ( v(k,j,i) +                   &
     2263                                                                      v(k,j+1,i) )                 &
     2264                                                         )**2 +                                    &
     2265                                                         ( 0.5_wp * ( w(k-1,j,i) +                 &
     2266                                                                      w(k,j,i) )                   &
     2267                                                         )**2                                      &
     2268                                                       ) *                                         &
     2269                                                   ( s(k,j,i) - leaf_surface_conc )
    22802270                   ENDDO
    22812271                ENDDO
    2282              ENDDO   
     2272             ENDDO
    22832273
    22842274
     
    22872277
    22882278             WRITE( message_string, * ) 'wrong component: ', component
    2289              CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
     2279             CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 )
    22902280
    22912281       END SELECT
     
    22942284
    22952285
    2296 !------------------------------------------------------------------------------!
     2286!--------------------------------------------------------------------------------------------------!
    22972287! Description:
    22982288! ------------
    2299 !> Calculation of the tendency terms, accounting for the effect of the plant
    2300 !> canopy on momentum and scalar quantities.
     2289!> Calculation of the tendency terms, accounting for the effect of the plant canopy on momentum and
     2290!> scalar quantities.
    23012291!>
    2302 !> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
    2303 !> (defined on scalar grid), as initialized in subroutine pcm_init.
    2304 !> The lad on the w-grid is vertically interpolated from the surrounding
    2305 !> lad_s. The upper boundary of the canopy is defined on the w-grid at
    2306 !> k = pch_index. Here, the lad is zero.
     2292!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0 (defined on scalar grid), as
     2293!> initialized in subroutine pcm_init.
     2294!> The lad on the w-grid is vertically interpolated from the surrounding lad_s. The upper boundary
     2295!> of the canopy is defined on the w-grid at k = pch_index. Here, the lad is zero.
    23072296!>
    2308 !> The canopy drag must be limited (previously accounted for by calculation of
    2309 !> a limiting canopy timestep for the determination of the maximum LES timestep
    2310 !> in subroutine timestep), since it is physically impossible that the canopy
    2311 !> drag alone can locally change the sign of a velocity component. This
    2312 !> limitation is realized by calculating preliminary tendencies and velocities.
    2313 !> It is subsequently checked if the preliminary new velocity has a different
    2314 !> sign than the current velocity. If so, the tendency is limited in a way that
    2315 !> the velocity can at maximum be reduced to zero by the canopy drag.
     2297!> The canopy drag must be limited (previously accounted for by calculation of a limiting canopy
     2298!> timestep for the determination of the maximum LES timestep in subroutine timestep), since it is
     2299!> physically impossible that the canopy drag alone can locally change the sign of a velocity
     2300!> component. This limitation is realized by calculating preliminary tendencies and velocities. It
     2301!> is subsequently checked if the preliminary new velocity has a different sign than the current
     2302!> velocity. If so, the tendency is limited in a way that the velocity can at maximum be reduced to
     2303!> zero by the canopy drag.
    23162304!>
    23172305!>
    23182306!> Call for grid point i,j
    2319 !------------------------------------------------------------------------------!
     2307!--------------------------------------------------------------------------------------------------!
    23202308    SUBROUTINE pcm_tendency_ij( i, j, component )
    23212309
     
    23352323       REAL(wp) ::  lad_local !< local lad value
    23362324       REAL(wp) ::  pre_tend  !< preliminary tendency
    2337        REAL(wp) ::  pre_u     !< preliminary u-value 
    2338        REAL(wp) ::  pre_v     !< preliminary v-value 
    2339        REAL(wp) ::  pre_w     !< preliminary w-value 
     2325       REAL(wp) ::  pre_u     !< preliminary u-value
     2326       REAL(wp) ::  pre_v     !< preliminary v-value
     2327       REAL(wp) ::  pre_w     !< preliminary w-value
    23402328
    23412329
     
    23492337          CASE ( 1 )
    23502338!
    2351 !--          Set control flags indicating east- and westward-orientated
    2352 !--          building edges. Note, building_egde_w is set from the perspective
    2353 !--          of the potential rooftop grid point, while building_edge_e is
    2354 !--          is set from the perspective of the non-building grid point.
    2355              building_edge_w = ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )  .AND. &
    2356                          .NOT. ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )
    2357              building_edge_e = ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )  .AND. &
    2358                          .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
     2339!--          Set control flags indicating east- and westward-orientated building edges. Note,
     2340!--          building_egde_w is set from the perspective of the potential rooftop grid point, while
     2341!--          building_edge_e is set from the perspective of the non-building grid point.
     2342             building_edge_w =       ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )  .AND.         &
     2343                               .NOT. ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )
     2344             building_edge_e =       ANY( BTEST( wall_flags_total_0(:,j,i-1), 6 ) )  .AND.         &
     2345                               .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
    23592346!
    23602347!--          Determine topography-top index on u-grid
     
    23642351
    23652352!
    2366 !--             In order to create sharp boundaries of the plant canopy,
    2367 !--             the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
    2368 !--             rather than being interpolated from the surrounding lad_s,
    2369 !--             because this would yield smaller lad at the canopy boundaries
     2353!--             In order to create sharp boundaries of the plant canopy, the lad on the u-grid at
     2354!--             index (k,j,i) is equal to lad_s(k,j,i), rather than being interpolated from the
     2355!--             surrounding lad_s, because this would yield smaller lad at the canopy boundaries
    23702356!--             than inside of the canopy.
    2371 !--             For the same reason, the lad at the rightmost(i+1)canopy
    2372 !--             boundary on the u-grid equals lad_s(k,j,i), which is considered
    2373 !--             in the next if-statement. Note, at left-sided building edges
    2374 !--             this is not applied, here the LAD is equals the LAD at grid
    2375 !--             point (k,j,i), in order to avoid that LAD is mistakenly mapped
    2376 !--             on top of a roof where (usually) is no LAD is defined.
     2357!--             For the same reason, the lad at the rightmost(i+1)canopy boundary on the u-grid
     2358!--             equals lad_s(k,j,i), which is considered in the next if-statement. Note, at
     2359!--             left-sided building edges this is not applied, here the LAD is equals the LAD at
     2360!--             grid point (k,j,i), in order to avoid that LAD is mistakenly mapped on top of a roof
     2361!--             where (usually) is no LAD is defined.
    23772362!--             The same is also valid for bad_s.
    23782363                lad_local = lad_s(kk,j,i)
     
    23842369                     .NOT. building_edge_w )  bad_local = bad_s(kk,j,i-1)
    23852370!
    2386 !--             In order to avoid that LAD is mistakenly considered at right-
    2387 !--             sided building edges (here the topography-top index for the
    2388 !--             u-component at index j,i is still on the building while the
    2389 !--             topography top for the scalar isn't), LAD is taken from grid
     2371!--             In order to avoid that LAD is mistakenly considered at right-sided building edges
     2372!--             (here the topography-top index for the u-component at index j,i is still on the
     2373!--             building while the topography top for the scalar isn't), LAD is taken from grid
    23902374!--             point (j,i-1). The same is also valid for bad_s.
    23912375                IF ( lad_local > 0.0_wp  .AND.  lad_s(kk,j,i-1) == 0.0_wp  .AND.                   &
     
    23982382!
    23992383!--             Calculate preliminary value (pre_tend) of the tendency
    2400                 pre_tend = - canopy_drag_coeff *                               &
    2401                              ( lad_local + bad_local ) *                       &
    2402                              SQRT( u(k,j,i)**2 +                               &
    2403                                    ( 0.25_wp * ( v(k,j,i-1)  +                 &
    2404                                                  v(k,j,i)    +                 &
    2405                                                  v(k,j+1,i)  +                 &
    2406                                                  v(k,j+1,i-1) )                &
    2407                                    )**2 +                                      &
    2408                                    ( 0.25_wp * ( w(k-1,j,i-1) +                &
    2409                                                  w(k-1,j,i)   +                &
    2410                                                  w(k,j,i-1)   +                &
    2411                                                  w(k,j,i) )                    &
    2412                                    )**2                                        &
    2413                                  ) *                                           &
     2384                pre_tend = - canopy_drag_coeff *                                                   &
     2385                             ( lad_local + bad_local ) *                                           &
     2386                             SQRT( u(k,j,i)**2 +                                                   &
     2387                                   ( 0.25_wp * ( v(k,j,i-1)  +                                     &
     2388                                                 v(k,j,i)    +                                     &
     2389                                                 v(k,j+1,i)  +                                     &
     2390                                                 v(k,j+1,i-1) )                                    &
     2391                                   )**2 +                                                          &
     2392                                   ( 0.25_wp * ( w(k-1,j,i-1) +                                    &
     2393                                                 w(k-1,j,i)   +                                    &
     2394                                                 w(k,j,i-1)   +                                    &
     2395                                                 w(k,j,i) )                                        &
     2396                                   )**2                                                            &
     2397                                 ) *                                                               &
    24142398                             u(k,j,i)
    24152399
     
    24182402                pre_u = u(k,j,i) + dt_3d * pre_tend
    24192403!
    2420 !--             Compare sign of old velocity and new preliminary velocity,
    2421 !--             and in case the signs are different, limit the tendency
    2422                 IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
     2404!--             Compare sign of old velocity and new preliminary velocity, and in case the signs are
     2405!--             different, limit the tendency.
     2406                IF ( SIGN( pre_u,u(k,j,i) ) /= pre_u )  THEN
    24232407                   pre_tend = - u(k,j,i) * ddt_3d
    24242408                ELSE
     
    24352419          CASE ( 2 )
    24362420!
    2437 !--          Set control flags indicating north- and southward-orientated
    2438 !--          building edges. Note, building_egde_s is set from the perspective
    2439 !--          of the potential rooftop grid point, while building_edge_n is
    2440 !--          is set from the perspective of the non-building grid point.
    2441              building_edge_s = ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )  .AND. &
    2442                          .NOT. ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )
    2443              building_edge_n = ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )  .AND. &
    2444                          .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
     2421!--          Set control flags indicating north- and southward-orientated building edges. Note,
     2422!--          building_egde_s is set from the perspective of the potential rooftop grid point, while
     2423!--          building_edge_n is set from the perspective of the non-building grid point.
     2424             building_edge_s =       ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )  .AND.         &
     2425                               .NOT. ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )
     2426             building_edge_n =       ANY( BTEST( wall_flags_total_0(:,j-1,i), 6 ) )  .AND.         &
     2427                               .NOT. ANY( BTEST( wall_flags_total_0(:,j,i),   6 ) )
    24452428!
    24462429!--          Determine topography-top index on v-grid
     
    24492432                kk = k - topo_top_ind(j,i,2)  !- lad arrays are defined flat
    24502433!
    2451 !--             In order to create sharp boundaries of the plant canopy,
    2452 !--             the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
    2453 !--             rather than being interpolated from the surrounding lad_s,
    2454 !--             because this would yield smaller lad at the canopy boundaries
     2434!--             In order to create sharp boundaries of the plant canopy, the lad on the v-grid at
     2435!--             index (k,j,i) is equal to lad_s(k,j,i), rather than being interpolated from the
     2436!--             surrounding lad_s, because this would yield smaller lad at the canopy boundaries
    24552437!--             than inside of the canopy.
    2456 !--             For the same reason, the lad at the northmost(j+1)canopy
    2457 !--             boundary on the v-grid equals lad_s(k,j,i), which is considered
    2458 !--             in the next if-statement. Note, at left-sided building edges
    2459 !--             this is not applied, here the LAD is equals the LAD at grid
    2460 !--             point (k,j,i), in order to avoid that LAD is mistakenly mapped
    2461 !--             on top of a roof where (usually) is no LAD is defined.
     2438!--             For the same reason, the lad at the northmost (j+1) canopy boundary on the v-grid
     2439!--             equals lad_s(k,j,i), which is considered in the next if-statement. Note, at
     2440!--             left-sided building edges this is not applied, here the LAD is equals the LAD at
     2441!--             grid point (k,j,i), in order to avoid that LAD is mistakenly mapped on top of a roof
     2442!--             where (usually) is no LAD is defined.
    24622443!--             The same is also valid for bad_s.
    24632444                lad_local = lad_s(kk,j,i)
     
    24692450                     .NOT. building_edge_s )  bad_local = bad_s(kk,j-1,i)
    24702451!
    2471 !--             In order to avoid that LAD is mistakenly considered at right-
    2472 !--             sided building edges (here the topography-top index for the
    2473 !--             u-component at index j,i is still on the building while the
    2474 !--             topography top for the scalar isn't), LAD is taken from grid
     2452!--             In order to avoid that LAD is mistakenly considered at right-sided building edges
     2453!--             (here the topography-top index for the u-component at index j,i is still on the
     2454!--             building while the topography top for the scalar isn't), LAD is taken from grid
    24752455!--             point (j,i-1). The same is also valid for bad_s.
    24762456                IF ( lad_local > 0.0_wp  .AND.  lad_s(kk,j-1,i) == 0.0_wp  .AND.                   &
     
    24832463!
    24842464!--             Calculate preliminary value (pre_tend) of the tendency
    2485                 pre_tend = - canopy_drag_coeff *                               &
    2486                              ( lad_local + bad_local ) *                       &
    2487                              SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                &
    2488                                                  u(k,j-1,i+1) +                &
    2489                                                  u(k,j,i)     +                &
    2490                                                  u(k,j,i+1) )                  &
    2491                                    )**2 +                                      &
    2492                                    v(k,j,i)**2 +                               &
    2493                                    ( 0.25_wp * ( w(k-1,j-1,i) +                &
    2494                                                  w(k-1,j,i)   +                &
    2495                                                  w(k,j-1,i)   +                &
    2496                                                  w(k,j,i) )                    &
    2497                                    )**2                                        &
    2498                                  ) *                                           &
     2465                pre_tend = - canopy_drag_coeff *                                                   &
     2466                             ( lad_local + bad_local ) *                                           &
     2467                             SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                                    &
     2468                                                 u(k,j-1,i+1) +                                    &
     2469                                                 u(k,j,i)     +                                    &
     2470                                                 u(k,j,i+1) )                                      &
     2471                                   )**2 +                                                          &
     2472                                   v(k,j,i)**2 +                                                   &
     2473                                   ( 0.25_wp * ( w(k-1,j-1,i) +                                    &
     2474                                                 w(k-1,j,i)   +                                    &
     2475                                                 w(k,j-1,i)   +                                    &
     2476                                                 w(k,j,i) )                                        &
     2477                                   )**2                                                            &
     2478                                 ) *                                                               &
    24992479                             v(k,j,i)
    25002480
     
    25052485!--             Compare sign of old velocity and new preliminary velocity,
    25062486!--             and in case the signs are different, limit the tendency
    2507                 IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
     2487                IF ( SIGN( pre_v,v(k,j,i) ) /= pre_v )  THEN
    25082488                   pre_tend = - v(k,j,i) * ddt_3d
    25092489                ELSE
     
    25292509!
    25302510!--             Calculate preliminary value (pre_tend) of the tendency
    2531                 pre_tend = - canopy_drag_coeff *                               &
    2532                              ( 0.5_wp *                                        &
    2533                                 ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) +          &
    2534                                0.5_wp *                                        &
    2535                                 ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) *        &
    2536                              SQRT( ( 0.25_wp * ( u(k,j,i)    +                 &
    2537                                                  u(k,j,i+1)  +                 &
    2538                                                  u(k+1,j,i)  +                 &
    2539                                                  u(k+1,j,i+1) )                &
    2540                                    )**2 +                                      &
    2541                                    ( 0.25_wp * ( v(k,j,i)    +                 &
    2542                                                  v(k,j+1,i)  +                 &
    2543                                                  v(k+1,j,i)  +                 &
    2544                                                  v(k+1,j+1,i) )                &
    2545                                    )**2 +                                      &
    2546                                    w(k,j,i)**2                                 &
    2547                                  ) *                                           &
     2511                pre_tend = - canopy_drag_coeff *                                                   &
     2512                             ( 0.5_wp * ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) +                      &
     2513                               0.5_wp * ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) *                    &
     2514                             SQRT( ( 0.25_wp * ( u(k,j,i)    +                                     &
     2515                                                 u(k,j,i+1)  +                                     &
     2516                                                 u(k+1,j,i)  +                                     &
     2517                                                 u(k+1,j,i+1) )                                    &
     2518                                   )**2 +                                                          &
     2519                                   ( 0.25_wp * ( v(k,j,i)    +                                     &
     2520                                                 v(k,j+1,i)  +                                     &
     2521                                                 v(k+1,j,i)  +                                     &
     2522                                                 v(k+1,j+1,i) )                                    &
     2523                                   )**2 +                                                          &
     2524                                   w(k,j,i)**2                                                     &
     2525                                 ) *                                                               &
    25482526                             w(k,j,i)
    25492527!
     
    25512529                pre_w = w(k,j,i) + dt_3d * pre_tend
    25522530!
    2553 !--             Compare sign of old velocity and new preliminary velocity,
    2554 !--             and in case the signs are different, limit the tendency
    2555                 IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
     2531!--             Compare sign of old velocity and new preliminary velocity, and in case the signs are
     2532!--             different, limit the tendency.
     2533                IF ( SIGN( pre_w,w(k,j,i) ) /= pre_w )  THEN
    25562534                   pre_tend = - w(k,j,i) * ddt_3d
    25572535                ELSE
     
    25682546!
    25692547!--          Determine topography-top index on scalar grid
    2570              IF ( humidity ) THEN
     2548             IF ( humidity )  THEN
    25712549                DO  k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
    25722550                   kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
    2573                    tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i) -       &
    2574                                  pcm_latent_rate(kk,j,i)
     2551                   tend(k,j,i) = tend(k,j,i) + pcm_heating_rate(kk,j,i) - pcm_latent_rate(kk,j,i)
    25752552                ENDDO
    25762553             ELSE
     
    25882565             DO  k = topo_top_ind(j,i,0) + 1, topo_top_ind(j,i,0) + pch_index_ji(j,i)
    25892566                kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
    2590                 IF ( .NOT. plant_canopy_transpiration ) THEN
     2567                IF ( .NOT. plant_canopy_transpiration )  THEN
    25912568                   ! pcm_transpiration_rate is calculated in radiation model
    25922569                   ! in case of plant_canopy_transpiration = .T.
    25932570                   ! to include also the dependecy to the radiation
    25942571                   ! in the plant canopy box
    2595                    pcm_transpiration_rate(kk,j,i) = - leaf_scalar_exch_coeff   &
    2596                                     * lad_s(kk,j,i) *                          &
    2597                                     SQRT( ( 0.5_wp * ( u(k,j,i) +              &
    2598                                                        u(k,j,i+1) )            &
    2599                                           )**2  +                              &
    2600                                           ( 0.5_wp * ( v(k,j,i) +              &
    2601                                                        v(k,j+1,i) )            &
    2602                                           )**2 +                               &
    2603                                           ( 0.5_wp * ( w(k-1,j,i) +            &
    2604                                                        w(k,j,i) )              &
    2605                                           )**2                                 &
    2606                                         ) *                                    &
    2607                                     ( q(k,j,i) - leaf_surface_conc )
     2572                   pcm_transpiration_rate(kk,j,i) = - leaf_scalar_exch_coeff                       &
     2573                                                    * lad_s(kk,j,i) *                              &
     2574                                                    SQRT( ( 0.5_wp * ( u(k,j,i) +                  &
     2575                                                                       u(k,j,i+1) )                &
     2576                                                          )**2  +                                  &
     2577                                                          ( 0.5_wp * ( v(k,j,i) +                  &
     2578                                                                       v(k,j+1,i) )                &
     2579                                                          )**2 +                                   &
     2580                                                          ( 0.5_wp * ( w(k-1,j,i) +                &
     2581                                                                       w(k,j,i) )                  &
     2582                                                          )**2                                     &
     2583                                                        ) *                                        &
     2584                                                    ( q(k,j,i) - leaf_surface_conc )
    26082585                ENDIF
    26092586
    26102587                tend(k,j,i) = tend(k,j,i) + pcm_transpiration_rate(kk,j,i)
    26112588
    2612              ENDDO   
     2589             ENDDO
    26132590
    26142591!
     
    26202597
    26212598                kk = k - topo_top_ind(j,i,0)
    2622                 tend(k,j,i) = tend(k,j,i) -                                    &
    2623                                  2.0_wp * canopy_drag_coeff *                  &
    2624                                  ( lad_s(kk,j,i) + bad_s(kk,j,i) ) *           &
    2625                                  SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
    2626                                                     u(k,j,i+1) )               &
    2627                                        )**2 +                                  & 
    2628                                        ( 0.5_wp * ( v(k,j,i) +                 &
    2629                                                     v(k,j+1,i) )               &
    2630                                        )**2 +                                  &
    2631                                        ( 0.5_wp * ( w(k,j,i) +                 &
    2632                                                     w(k+1,j,i) )               &
    2633                                        )**2                                    &
    2634                                      ) *                                       &
    2635                                  e(k,j,i)
     2599                tend(k,j,i) = tend(k,j,i) -  2.0_wp * canopy_drag_coeff *                          &
     2600                                             ( lad_s(kk,j,i) + bad_s(kk,j,i) ) *                   &
     2601                                             SQRT( ( 0.5_wp * ( u(k,j,i) +                         &
     2602                                                                u(k,j,i+1) )                       &
     2603                                                   )**2 +                                          &
     2604                                                   ( 0.5_wp * ( v(k,j,i) +                         &
     2605                                                                v(k,j+1,i) )                       &
     2606                                                   )**2 +                                          &
     2607                                                   ( 0.5_wp * ( w(k,j,i) +                         &
     2608                                                                w(k+1,j,i) )                       &
     2609                                                   )**2                                            &
     2610                                                 ) *                                               &
     2611                                             e(k,j,i)
    26362612             ENDDO
    26372613!
    2638 !--       scalar concentration 
     2614!--       scalar concentration
    26392615          CASE ( 7 )
    26402616!
     
    26432619
    26442620                kk = k - topo_top_ind(j,i,0)
    2645                 tend(k,j,i) = tend(k,j,i) -                                    &
    2646                                  leaf_scalar_exch_coeff *                      &
    2647                                  lad_s(kk,j,i) *                               &
    2648                                  SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
    2649                                                     u(k,j,i+1) )               &
    2650                                        )**2  +                                 &
    2651                                        ( 0.5_wp * ( v(k,j,i) +                 &
    2652                                                     v(k,j+1,i) )               &
    2653                                        )**2 +                                  &
    2654                                        ( 0.5_wp * ( w(k-1,j,i) +               &
    2655                                                     w(k,j,i) )                 &
    2656                                        )**2                                    &
    2657                                      ) *                                       &
    2658                                  ( s(k,j,i) - leaf_surface_conc )
    2659              ENDDO               
     2621                tend(k,j,i) = tend(k,j,i) -  leaf_scalar_exch_coeff *                              &
     2622                                             lad_s(kk,j,i) *                                       &
     2623                                             SQRT( ( 0.5_wp * ( u(k,j,i) +                         &
     2624                                                                u(k,j,i+1) )                       &
     2625                                                   )**2  +                                         &
     2626                                                   ( 0.5_wp * ( v(k,j,i) +                         &
     2627                                                                v(k,j+1,i) )                       &
     2628                                                   )**2 +                                          &
     2629                                                   ( 0.5_wp * ( w(k-1,j,i) +                       &
     2630                                                                w(k,j,i) )                         &
     2631                                                   )**2                                            &
     2632                                                 ) *                                               &
     2633                                             ( s(k,j,i) - leaf_surface_conc )
     2634             ENDDO
    26602635
    26612636       CASE DEFAULT
    26622637
    26632638          WRITE( message_string, * ) 'wrong component: ', component
    2664           CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
     2639          CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 )
    26652640
    26662641       END SELECT
     
    26682643    END SUBROUTINE pcm_tendency_ij
    26692644
    2670 !------------------------------------------------------------------------------!
     2645!--------------------------------------------------------------------------------------------------!
    26712646! Description:
    26722647! ------------
    26732648!> Subroutine writes global restart data
    2674 !------------------------------------------------------------------------------!
     2649!--------------------------------------------------------------------------------------------------!
    26752650    SUBROUTINE pcm_wrd_global
    26762651
     
    26782653
    26792654          CALL wrd_write_string( 'pch_index' )
    2680           WRITE ( 14 )  pch_index
    2681 
    2682        ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
     2655          WRITE( 14 )  pch_index
     2656
     2657       ELSE IF ( restart_data_format_output(1:3) == 'mpi' )  THEN
    26832658
    26842659          CALL wrd_mpi_io( 'pch_index', pch_index )
     
    26882663    END SUBROUTINE pcm_wrd_global
    26892664
    2690 !------------------------------------------------------------------------------!
     2665!--------------------------------------------------------------------------------------------------!
    26912666! Description:
    26922667! ------------
    26932668!> Subroutine writes local (subdomain) restart data
    2694 !------------------------------------------------------------------------------!
     2669!--------------------------------------------------------------------------------------------------!
    26952670    SUBROUTINE pcm_wrd_local
    26962671
     
    26982673                                                           !< non-standard vertical index bounds
    26992674
     2675
    27002676       IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
    27012677
    27022678          IF ( ALLOCATED( pcm_heatrate_av ) )  THEN
    27032679             CALL wrd_write_string( 'pcm_heatrate_av' )
    2704              WRITE ( 14 )  pcm_heatrate_av
     2680             WRITE( 14 )  pcm_heatrate_av
    27052681          ENDIF
    27062682
    27072683          IF ( ALLOCATED( pcm_latentrate_av ) )  THEN
    27082684             CALL wrd_write_string( 'pcm_latentrate_av' )
    2709              WRITE ( 14 )  pcm_latentrate_av
     2685             WRITE( 14 )  pcm_latentrate_av
    27102686          ENDIF
    27112687
    27122688          IF ( ALLOCATED( pcm_transpirationrate_av ) )  THEN
    27132689             CALL wrd_write_string( 'pcm_transpirationrate_av' )
    2714              WRITE ( 14 )  pcm_transpirationrate_av
     2690             WRITE( 14 )  pcm_transpirationrate_av
    27152691          ENDIF
    27162692
    2717        ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
     2693       ELSE IF ( restart_data_format_output(1:3) == 'mpi' )  THEN
    27182694
    27192695!
Note: See TracChangeset for help on using the changeset viewer.