Home

Context Navigation

← Previous Changeset
Next Changeset →

Changeset 4624

Timestamp:

Jul 24, 2020 9:53:17 AM (6 months ago)

Author:

raasch

Message:

file re-formatted to follow the PALM coding standard

File:

: 1 edited

palm/trunk/SOURCE/radiation_model_mod.f90 (modified) (235 diffs)

Legend:

: Unmodified
: Added
: Removed

palm/trunk/SOURCE/radiation_model_mod.f90

-                      r4587
+                      r4624
 !> @file radiation_model_mod.f90
 !------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! This file is part of the PALM model system.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 2015-2020 Institute of Computer Science of the
+!                     Czech Academy of Sciences, Prague
+! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
+! Public License as published by the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
+! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+! Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with PALM. If not, see
+! <http://www.gnu.org/licenses/>.
+!
+! Copyright 2015-2020 Institute of Computer Science of the Czech Academy of Sciences, Prague
 ! Copyright 2015-2020 Czech Technical University in Prague
 ! Copyright 1997-2020 Leibniz Universitaet Hannover
+!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
+!
+!
 ! Current revisions:
 ! ------------------
+! -----------------
+!
+!
 …
 ! -----------------
 ! $Id$
+! File re-formatted to follow the PALM coding standard
+!
+! 4587 2020-07-06 08:53:45Z pavelkrc
 ! RTM version 3.1 (see also previous commits):
 ! - rotation_angle supported
 ! - plant canopy box count minimized
 ! - multiple enhancements and bugfixes
+!
+! - Rotation_angle supported
+! - Plant canopy box count minimized
+! - Multiple enhancements and bugfixes
+!
 ! 4584 2020-06-29 13:16:14Z pavelkrc
 ! Consider only boxes with LAD>0 as plant canopy (credit: S. Schubert)
+!
+!
 ! 4576 2020-06-24 17:58:55Z pavelkrc
 ! Allow the use of rotation_angle in RTM
+!
+!
 ! 4574 2020-06-24 16:33:32Z pavelkrc
 ! - Restructure code in radiation_check_data_output
 ! - Move calculation of MPI global array offsets to a subroutine
+!
+!
 ! 4571 2020-06-24 08:59:06Z sebschub
 ! Bugfix in vertical lad_s coordinate
+!
+!
 ! 4558 2020-06-10 16:27:30Z moh.hefny
 ! Bugfix: - reset RTM output average values after each averaging timestep to zero
 !         - correct calculation of rtm_rad_net_av
+!
+! Bugfix: - Reset RTM output average values after each averaging timestep to zero
+!         - Correct calculation of rtm_rad_net_av
+!
 ! 4555 2020-06-05 21:52:00Z moh.hefny
 ! Bugfix in averaging PC and MRT related quantities
+!
+!
 ! 4552 2020-06-02 20:33:29Z moh.hefny
 ! Bugfix in IF statement in the emissivity coupling parameter for radiation-RTM
+!
+!
 ! 4535 2020-05-15 12:07:23Z raasch
 ! bugfix for restart data format query
+!
+! Bugfix for restart data format query
+!
 ! 4531 2020-05-13 09:52:22Z moh.hefny
 ! Bugfix in gather flux pabs_pc_lwdif in non_parallel case
+!
 ! 4529 2020-05-12 09:14:57Z moh.hefny
 ! - added the following new features to the coupling of RTM-radiation model:
 !   1) considering the vegetation interaction with LW in the coupling
 !   2) considering PC emissivity in calculating the effective emissivity
 !   3) new algorithm for claculating the coupling parameters so that each term
 !      is calculated within its original line and not at the end.
 ! - minor formatting and comments changes
+!
+! - Added the following new features to the coupling of RTM-radiation model:
+!   1) Considering the vegetation interaction with LW in the coupling
+!   2) Considering PC emissivity in calculating the effective emissivity
+!   3) New algorithm for claculating the coupling parameters so that each term is calculated within
+!      its original line and not at the end.
+! - Minor formatting and comments changes
+!
 ! 4517 2020-05-03 14:29:30Z raasch
 ! added restart with MPI-IO for reading local arrays
+!
+! Added restart with MPI-IO for reading local arrays
+!
 ! 4495 2020-04-13 20:11:20Z raasch
 ! restart data handling with MPI-IO added
+!
+! Restart data handling with MPI-IO added
+!
 ! 4493 2020-04-10 09:49:43Z pavelkrc
 ! Avoid unstable direct normal radiation near horizon
+!
+!
 ! 4481 2020-03-31 18:55:54Z maronga
 ! use statement for exchange horiz added
+!
+! Use statement for exchange horiz added
+!
 ! 4452 2020-03-10 20:15:32Z suehring
 ! Bugfix in calc_albedo
+!
 ! 4442 2020-03-04 19:21:13Z suehring
 ! - Change order of dimension in surface arrays %frac, %emissivity and %albedo
 !   to allow for better vectorization in the radiation interactions.
+! - Change order of dimension in surface arrays %frac, %emissivity and %albedo to allow for better
+!   vectorization in the radiation interactions.
 ! - Minor formatting issues
+!
 ! 4441 2020-03-04 19:20:35Z suehring
 ! bugfixes: cpp-directives for serial mode moved, small changes to get serial mode compiled
+! Bugfixes: cpp-directives for serial mode moved, small changes to get serial mode compiled
+!
 ! 4400 2020-02-10 20:32:41Z suehring
 …
+!
 ! 4360 2020-01-07 11:25:50Z suehring
 ! Renamed pc_heating_rate, pc_transpiration_rate, pc_transpiration_rate to
 ! pcm_heating_rate, pcm_latent_rate, pcm_transpiration_rate
+! Renamed pc_heating_rate, pc_transpiration_rate, pc_transpiration_rate to pcm_heating_rate,
+! pcm_latent_rate, pcm_transpiration_rate
+!
 ! 4340 2019-12-16 08:17:03Z Giersch
 ! Albedo indices for building_surface_pars are now declared as parameters to
 ! prevent an error if the gfortran compiler with -Werror=unused-value is used
+! Albedo indices for building_surface_pars are now declared as parameters to prevent an error if the
+! gfortran compiler with -Werror=unused-value is used
+!
 ! 4291 2019-11-11 12:36:54Z moh.hefny
 ! Enabled RTM in case of biometeorology even if there is no vertical
 ! surfaces or 3D vegetation in the domain
+! Enabled RTM in case of biometeorology even if there is no vertical surfaces or 3D vegetation in
+! the domain
+!
 ! 4286 2019-10-30 16:01:14Z resler
 …
+!
 ! 4227 2019-09-10 18:04:34Z gronemeier
 ! implement new palm_date_time_mod
+! Implement new palm_date_time_mod
+!
 ! 4226 2019-09-10 17:03:24Z suehring
 …
 ! - Revise steering of splitting diffuse and direct radiation
 ! - Bugfixes in checks
+! - Optimize mapping of radiation components onto 2D arrays, avoid unnecessary
+!   operations
+! - Optimize mapping of radiation components onto 2D arrays, avoid unnecessary operations
+!
 ! 4208 2019-09-02 09:01:07Z suehring
+! Bugfix in accessing albedo_pars in the clear-sky branch
+! (merge from branch resler)
+! Bugfix in accessing albedo_pars in the clear-sky branch (merge from branch resler)
+!
 ! 4198 2019-08-29 15:17:48Z gronemeier
 …
+!
 ! 4190 2019-08-27 15:42:37Z suehring
+! Implement external radiation forcing also for level-of-detail = 2
+! (horizontally 2D radiation)
+! Implement external radiation forcing also for level-of-detail = 2 (horizontally 2D radiation)
+!
 ! 4188 2019-08-26 14:15:47Z suehring
 …
+!
 ! 4187 2019-08-26 12:43:15Z suehring
+! - Take external radiation from root domain dynamic input if not provided for
+!   each nested domain
+! - Take external radiation from root domain dynamic input if not provided for each nested domain
 ! - Combine MPI_ALLREDUCE calls to reduce mpi overhead
+!
 …
+!
 ! 4089 2019-07-11 14:30:27Z suehring
 ! - Correct level 2 initialization of spectral albedos in rrtmg branch, long- and
 !   shortwave albedos were mixed-up.
 ! - Change order of albedo_pars so that it is now consistent with the defined
 !   order of albedo_pars in PIDS
+! - Correct level 2 initialization of spectral albedos in rrtmg branch, long- and shortwave albedos
+!   were mixed-up.
+! - Change order of albedo_pars so that it is now consistent with the defined order of albedo_pars
+!   in PIDS
+!
 ! 4069 2019-07-01 14:05:51Z Giersch
 ! Masked output running index mid has been introduced as a local variable to
 ! avoid runtime error (Loop variable has been modified) in time_integration
+! Masked output running index mid has been introduced as a local variable to avoid runtime error
+! (Loop variable has been modified) in time_integration
+!
 ! 4067 2019-07-01 13:29:25Z suehring
 …
+!
 ! 4008 2019-05-30 09:50:11Z moh.hefny
+! Bugfix in check variable when a variable's string is less than 3
+! characters is processed. All variables now are checked if they
+! belong to radiation
+! Bugfix in check variable when a variable's string is less than 3 characters is processed. All
+! variables now are checked if they belong to radiation
+!
 ! 3992 2019-05-22 16:49:38Z suehring
 ! Bugfix in rrtmg radiation branch in a nested run when the lowest prognistic
 ! grid points in a child domain are all inside topography
+! Bugfix in rrtmg radiation branch in a nested run when the lowest prognistic grid points in a child
+! domain are all inside topography
+!
 ! 3987 2019-05-22 09:52:13Z kanani
 …
+!
 ! 3885 2019-04-11 11:29:34Z kanani
 ! Changes related to global restructuring of location messages and introduction
 ! of additional debug messages
+! Changes related to global restructuring of location messages and introduction of additional debug
+! messages
+!
 ! 3881 2019-04-10 09:31:22Z suehring
+! Output of albedo and emissivity moved from USM, bugfixes in initialization
+! of albedo
+! Output of albedo and emissivity moved from USM, bugfixes in initialization of albedo
+!
 ! 3861 2019-04-04 06:27:41Z maronga
 …
+!
 ! 3846 2019-04-01 13:55:30Z suehring
 ! unused variable removed
+! Unused variable removed
+!
 ! 3814 2019-03-26 08:40:31Z pavelkrc
 …
+!
 ! 3771 2019-02-28 12:19:33Z raasch
 ! rrtmg preprocessor for directives moved/added, save attribute added to temporary
 ! pointers to avoid compiler warnings about outlived pointer targets,
 ! statement added to avoid compiler warning about unused variable
+! rrtmg preprocessor for directives moved/added, save attribute added to temporary pointers to avoid
+! compiler warnings about outlived pointer targets, statement added to avoid compiler warning about
+! unused variable
+!
 ! 3769 2019-02-28 10:16:49Z moh.hefny
 ! removed unused variables and subroutine radiation_radflux_gridbox
+! Removed unused variables and subroutine radiation_radflux_gridbox
+!
 ! 3767 2019-02-27 08:18:02Z raasch
 ! unused variable for file index removed from rrd-subroutines parameter list
+! Unused variable for file index removed from rrd-subroutines parameter list
+!
 ! 3760 2019-02-21 18:47:35Z moh.hefny
 ! Bugfix: initialized simulated_time before calculating solar position
 ! to enable restart option with reading in SVF from file(s).
+! Bugfix: initialized simulated_time before calculating solar position to enable restart option with
+! reading in SVF from file(s).
+!
 ! 3754 2019-02-19 17:02:26Z kanani
 ! (resler, pavelkrc)
+! Bugfixes: add further required MRT factors to read/write_svf,
+! fix for aggregating view factors to eliminate local noise in reflected
+! irradiance at mutually close surfaces (corners, presence of trees) in the
+! angular discretization scheme.
+! Bugfixes: add further required MRT factors to read/write_svf, fix for aggregating view factors to
+! eliminate local noise in reflected irradiance at mutually close surfaces (corners, presence of
+! trees) in the angular discretization scheme.
+!
 ! 3752 2019-02-19 09:37:22Z resler
 ! added read/write number of MRT factors to the respective routines
+! Added read/write number of MRT factors to the respective routines
+!
 ! 3705 2019-01-29 19:56:39Z suehring
 …
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! nopointer option removed
+! Nopointer option removed
+!
 ! 1496 2014-12-02 17:25:50Z maronga
 …
+!
+!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Radiation models and interfaces:
 !> constant, simple and RRTMG models, interface to external radiation model
+!> Constant, simple and RRTMG models, interface to external radiation model
 !> Radiative Transfer Model (RTM) version 3.0 for modelling of radiation
 !> interactions within urban canopy or other surface layer in complex terrain
+!> Interactions within urban canopy or other surface layer in complex terrain
 !> Integrations of RTM with other PALM-4U modules:
 !> integration with RRTMG, USM, LSM, PCM, BIO modules
+!> Integration with RRTMG, USM, LSM, PCM, BIO modules
 !>
 !> @todo move variable definitions used in radiation_init only to the subroutine
 !>       as they are no longer required after initialization.
+!> @todo Move variable definitions used in radiation_init only to the subroutine as they are no
+!>       longer required after initialization.
 !> @todo Output of full column vertical profiles used in RRTMG
 !> @todo Output of other rrtm arrays (such as volume mixing ratios)
 !> @todo Optimize radiation_tendency routines
 !>
 !> @note Many variables have a leading dummy dimension (0:0) in order to
 !>       match the assume-size shape expected by the RRTMG model.
 !------------------------------------------------------------------------------!
+!> @note Many variables have a leading dummy dimension (0:0) in order to match the assume-size shape
+!>       expected by the RRTMG model.
+!--------------------------------------------------------------------------------------------------!
  MODULE radiation_model_mod
+    USE arrays_3d,                                                             &
+        ONLY:  dzw, hyp, nc, pt, p, q, ql, u, v, w, zu, zw, exner, d_exner
+    USE basic_constants_and_equations_mod,                                     &
+        ONLY:  c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, sigma_sb, &
+               barometric_formula
+    USE calc_mean_profile_mod,                                                 &
+    USE arrays_3d,                                                                                 &
+        ONLY:  dzw,                                                                                &
+               d_exner,                                                                            &
+               exner,                                                                              &
+               hyp,                                                                                &
+               nc,                                                                                 &
+               pt,                                                                                 &
+               p,                                                                                  &
+               q,                                                                                  &
+               ql,                                                                                 &
+               u,                                                                                  &
+               v,                                                                                  &
+               w,                                                                                  &
+               zu,                                                                                 &
+               zw
+    USE basic_constants_and_equations_mod,                                                         &
+        ONLY:  barometric_formula,                                                                 &
+               c_p,                                                                                &
+               g,                                                                                  &
+               lv_d_cp,                                                                            &
+               l_v,                                                                                &
+               pi,                                                                                 &
+               r_d,                                                                                &
+               rho_l,                                                                              &
+               sigma_sb,                                                                           &
+               solar_constant
+    USE calc_mean_profile_mod,                                                                     &
         ONLY:  calc_mean_profile
+    USE control_parameters,                                                    &
+        ONLY:  biometeorology, cloud_droplets, coupling_char,                  &
+               debug_output, debug_output_timestep, debug_string,              &
+               dt_3d,                                                          &
+               dz, dt_spinup, end_time,                                        &
+               humidity,                                                       &
+               initializing_actions, io_blocks, io_group,                      &
+               land_surface, large_scale_forcing,                              &
+               latitude, longitude, lsf_surf,                                  &
+               message_string, plant_canopy, pt_surface,                       &
+               rho_surface, simulated_time, spinup_time, surface_pressure,     &
+               read_svf, restart_data_format_output, write_svf,                &
+               time_since_reference_point, urban_surface, varnamelength
+    USE cpulog,                                                                &
+        ONLY:  cpu_log, log_point, log_point_s
+    USE grid_variables,                                                        &
+         ONLY:  ddx, ddy, dx, dy
+    USE indices,                                                               &
+        ONLY:  nnx, nny, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,   &
+               nzb, nzt, topo_top_ind
+    USE control_parameters,                                                                        &
+        ONLY:  biometeorology,                                                                     &
+               cloud_droplets,                                                                     &
+               coupling_char,                                                                      &
+               debug_output,                                                                       &
+               debug_output_timestep,                                                              &
+               debug_string,                                                                       &
+               dt_3d,                                                                              &
+               dz,                                                                                 &
+               dt_spinup,                                                                          &
+               end_time,                                                                           &
+               humidity,                                                                           &
+               initializing_actions,                                                               &
+               io_blocks,                                                                          &
+               io_group,                                                                           &
+               land_surface,                                                                       &
+               large_scale_forcing,                                                                &
+               latitude,                                                                           &
+               longitude,                                                                          &
+               lsf_surf,                                                                           &
+               message_string,                                                                     &
+               plant_canopy,                                                                       &
+               pt_surface,                                                                         &
+               read_svf,                                                                           &
+               restart_data_format_output,                                                         &
+               rho_surface,                                                                        &
+               simulated_time,                                                                     &
+               spinup_time,                                                                        &
+               surface_pressure,                                                                   &
+               time_since_reference_point,                                                         &
+               urban_surface,                                                                      &
+               varnamelength,                                                                      &
+               write_svf
+    USE cpulog,                                                                                    &
+        ONLY:  cpu_log,                                                                            &
+               log_point,                                                                          &
+               log_point_s
+    USE grid_variables,                                                                            &
+         ONLY:  ddx,                                                                               &
+                ddy,                                                                               &
+                dx,                                                                                &
+                dy
+    USE indices,                                                                                   &
+        ONLY:  nnx,                                                                                &
+               nny,                                                                                &
+               nx,                                                                                 &
+               nxl,                                                                                &
+               nxlg,                                                                               &
+               nxr,                                                                                &
+               nxrg,                                                                               &
+               ny,                                                                                 &
+               nyn,                                                                                &
+               nyng,                                                                               &
+               nys,                                                                                &
+               nysg,                                                                               &
+               nzb,                                                                                &
+               nzt,                                                                                &
+               topo_top_ind
     USE, INTRINSIC :: iso_c_binding
 …
     USE kinds
+    USE bulk_cloud_model_mod,                                                  &
+        ONLY:  bulk_cloud_model, microphysics_morrison, na_init, nc_const, sigma_gc
+    USE bulk_cloud_model_mod,                                                                      &
+        ONLY:  bulk_cloud_model,                                                                   &
+               microphysics_morrison,                                                              &
+               na_init,                                                                            &
+               nc_const,                                                                           &
+               sigma_gc
 #if defined ( __netcdf )
 …
 #endif
+    USE netcdf_data_input_mod,                                                 &
+        ONLY:  albedo_type_f,                                                  &
+               albedo_pars_f,                                                  &
+               building_type_f,                                                &
+               building_surface_pars_f,                                        &
+               pavement_type_f,                                                &
+               vegetation_type_f,                                              &
+               water_type_f,                                                   &
+               char_fill,                                                      &
+               char_lod,                                                       &
+               check_existence,                                                &
+               close_input_file,                                               &
+               get_attribute,                                                  &
+               get_dimension_length,                                           &
+               get_variable,                                                   &
+               inquire_num_variables,                                          &
+               inquire_variable_names,                                         &
+               input_file_dynamic,                                             &
+               input_pids_dynamic,                                             &
+               num_var_pids,                                                   &
+               pids_id,                                                        &
+               open_read_file,                                                 &
+               real_1d_3d,                                                     &
+               vars_pids
+    USE palm_date_time_mod,                                                    &
+        ONLY:  date_time_str_len, get_date_time,                               &
+               hours_per_day, seconds_per_hour
+    USE plant_canopy_model_mod,                                                &
+        ONLY:  lad_s,                                                          &
+               pcm_heating_rate,                                               &
+               pcm_transpiration_rate,                                         &
+               pcm_latent_rate,                                                &
+               plant_canopy_transpiration,                                     &
+               pcm_calc_transpiration_rate
+    USE netcdf_data_input_mod,                                                                     &
+        ONLY:  albedo_type_f,                                                                      &
+               albedo_pars_f,                                                                      &
+               building_type_f,                                                                    &
+               building_surface_pars_f,                                                            &
+               char_fill,                                                                          &
+               char_lod,                                                                           &
+               check_existence,                                                                    &
+               close_input_file,                                                                   &
+               get_attribute,                                                                      &
+               get_dimension_length,                                                               &
+               get_variable,                                                                       &
+               inquire_num_variables,                                                              &
+               inquire_variable_names,                                                             &
+               input_file_dynamic,                                                                 &
+               input_pids_dynamic,                                                                 &
+               num_var_pids,                                                                       &
+               open_read_file,                                                                     &
+               pavement_type_f,                                                                    &
+               pids_id,                                                                            &
+               real_1d_3d,                                                                         &
+               vars_pids,                                                                          &
+               vegetation_type_f,                                                                  &
+               water_type_f
+    USE palm_date_time_mod,                                                                        &
+        ONLY:  date_time_str_len,                                                                  &
+               get_date_time,                                                                      &
+               hours_per_day,                                                                      &
+               seconds_per_hour
+    USE plant_canopy_model_mod,                                                                    &
+        ONLY:  lad_s,                                                                              &
+               pcm_calc_transpiration_rate,                                                        &
+               pcm_heating_rate,                                                                   &
+               pcm_transpiration_rate,                                                             &
+               pcm_latent_rate,                                                                    &
+               plant_canopy_transpiration
     USE pegrid
 #if defined ( __rrtmg )
+    USE parrrsw,                                                               &
+        ONLY:  naerec, nbndsw
+    USE parrrtm,                                                               &
+    USE parrrsw,                                                                                   &
+        ONLY:  naerec,                                                                             &
+               nbndsw
+    USE parrrtm,                                                                                   &
         ONLY:  nbndlw
     USE rrtmg_lw_init,                                                         &
+    USE rrtmg_lw_init,                                                                             &
         ONLY:  rrtmg_lw_ini
     USE rrtmg_sw_init,                                                         &
+    USE rrtmg_sw_init,                                                                             &
         ONLY:  rrtmg_sw_ini
     USE rrtmg_lw_rad,                                                          &
+    USE rrtmg_lw_rad,                                                                              &
         ONLY:  rrtmg_lw
     USE rrtmg_sw_rad,                                                          &
+    USE rrtmg_sw_rad,                                                                              &
         ONLY:  rrtmg_sw
 #endif
     USE restart_data_mpi_io_mod,                                                                   &
+        ONLY:  rd_mpi_io_check_array, rrd_mpi_io, wrd_mpi_io
+    USE statistics,                                                            &
+        ONLY:  rd_mpi_io_check_array,                                                              &
+               rrd_mpi_io,                                                                         &
+               wrd_mpi_io
+    USE statistics,                                                                                &
         ONLY:  hom
+    USE surface_mod,                                                           &
+        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win,                       &
+               surf_lsm_h, surf_lsm_v, surf_type, surf_usm_h, surf_usm_v,      &
+    USE surface_mod,                                                                               &
+        ONLY:  ind_pav_green,                                                                      &
+               ind_veg_wall,                                                                       &
+               ind_wat_win,                                                                        &
+               surf_lsm_h,                                                                         &
+               surf_lsm_v,                                                                         &
+               surf_type,                                                                          &
+               surf_usm_h,                                                                         &
+               surf_usm_v,                                                                         &
                vertical_surfaces_exist
     IMPLICIT NONE
     CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
+    CHARACTER(10) ::  radiation_scheme = 'clear-sky'  ! 'constant', 'clear-sky', or 'rrtmg'
+!
 !-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
     CHARACTER(37), DIMENSION(0:33), PARAMETER :: albedo_type_name = (/      &
                                    'user defined                         ', & !  0
                                    'ocean                                ', & !  1
                                    'mixed farming, tall grassland        ', & !  2
                                    'tall/medium grassland                ', & !  3
                                    'evergreen shrubland                  ', & !  4
                                    'short grassland/meadow/shrubland     ', & !  5
                                    'evergreen needleleaf forest          ', & !  6
                                    'mixed deciduous evergreen forest     ', & !  7
                                    'deciduous forest                     ', & !  8
                                    'tropical evergreen broadleaved forest', & !  9
                                    'medium/tall grassland/woodland       ', & ! 10
                                    'desert, sandy                        ', & ! 11
                                    'desert, rocky                        ', & ! 12
                                    'tundra                               ', & ! 13
                                    'land ice                             ', & ! 14
                                    'sea ice                              ', & ! 15
                                    'snow                                 ', & ! 16
                                    'bare soil                            ', & ! 17
                                    'asphalt/concrete mix                 ', & ! 18
                                    'asphalt (asphalt concrete)           ', & ! 19
                                    'concrete (Portland concrete)         ', & ! 20
                                    'sett                                 ', & ! 21
                                    'paving stones                        ', & ! 22
                                    'cobblestone                          ', & ! 23
                                    'metal                                ', & ! 24
                                    'wood                                 ', & ! 25
                                    'gravel                               ', & ! 26
                                    'fine gravel                          ', & ! 27
                                    'pebblestone                          ', & ! 28
                                    'woodchips                            ', & ! 29
                                    'tartan (sports)                      ', & ! 30
                                    'artifical turf (sports)              ', & ! 31
                                    'clay (sports)                        ', & ! 32
                                    'building (dummy)                     '  & ! 33
+    CHARACTER(37), DIMENSION(0:33), PARAMETER :: albedo_type_name = (/                       &
+                                   'user defined                         ',                  &  !<  0
+                                   'ocean                                ',                  &  !<  1
+                                   'mixed farming, tall grassland        ',                  &  !<  2
+                                   'tall/medium grassland                ',                  &  !<  3
+                                   'evergreen shrubland                  ',                  &  !<  4
+                                   'short grassland/meadow/shrubland     ',                  &  !<  5
+                                   'evergreen needleleaf forest          ',                  &  !<  6
+                                   'mixed deciduous evergreen forest     ',                  &  !<  7
+                                   'deciduous forest                     ',                  &  !<  8
+                                   'tropical evergreen broadleaved forest',                  &  !<  9
+                                   'medium/tall grassland/woodland       ',                  &  !< 10
+                                   'desert, sandy                        ',                  &  !< 11
+                                   'desert, rocky                        ',                  &  !< 12
+                                   'tundra                               ',                  &  !< 13
+                                   'land ice                             ',                  &  !< 14
+                                   'sea ice                              ',                  &  !< 15
+                                   'snow                                 ',                  &  !< 16
+                                   'bare soil                            ',                  &  !< 17
+                                   'asphalt/concrete mix                 ',                  &  !< 18
+                                   'asphalt (asphalt concrete)           ',                  &  !< 19
+                                   'concrete (Portland concrete)         ',                  &  !< 20
+                                   'sett                                 ',                  &  !< 21
+                                   'paving stones                        ',                  &  !< 22
+                                   'cobblestone                          ',                  &  !< 23
+                                   'metal                                ',                  &  !< 24
+                                   'wood                                 ',                  &  !< 25
+                                   'gravel                               ',                  &  !< 26
+                                   'fine gravel                          ',                  &  !< 27
+                                   'pebblestone                          ',                  &  !< 28
+                                   'woodchips                            ',                  &  !< 29
+                                   'tartan (sports)                      ',                  &  !< 30
+                                   'artifical turf (sports)              ',                  &  !< 31
+                                   'clay (sports)                        ',                  &  !< 32
+                                   'building (dummy)                     '                   &  !< 33
                                                          /)
+!
+!-- Indices of radiation-related input attributes in building_surface_pars
+!-- (other are in urban_surface_mod)
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_b_wall                = 19 !< index for Broadband albedo of wall fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_l_wall                = 20 !< index for Longwave albedo of wall fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_s_wall                = 21 !< index for Shortwave albedo of wall fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_b_win                 = 22 !< index for Broadband albedo of window fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_l_win                 = 23 !< index for Longwave albedo of window fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_s_win                 = 24 !< index for Shortwave albedo of window fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_b_green               = 24 !< index for Broadband albedo of green fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_l_green               = 25 !< index for Longwave albedo of green fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_s_green               = 26 !< index for Shortwave albedo of green fraction
+    INTEGER(iwp) :: albedo_type  = 9999999_iwp, &     !< Albedo surface type
+                    dots_rad     = 0_iwp              !< starting index for timeseries output
+    LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
+                constant_albedo = .FALSE.,            & !< flag parameter indicating whether the albedo may change depending on zenith
+                force_radiation_call = .FALSE.,       & !< flag parameter for unscheduled radiation calls
+                lw_radiation = .TRUE.,                & !< flag parameter indicating whether longwave radiation shall be calculated
+                radiation = .FALSE.,                  & !< flag parameter indicating whether the radiation model is used
+                sun_up    = .TRUE.,                   & !< flag parameter indicating whether the sun is up or down
+                sw_radiation = .TRUE.,                & !< flag parameter indicating whether shortwave radiation shall be calculated
+                sun_direction = .FALSE.,              & !< flag parameter indicating whether solar direction shall be calculated
+                average_radiation = .FALSE.,          & !< flag to set the calculation of radiation averaging for the domain
+                radiation_interactions = .FALSE.,     & !< flag to activiate RTM (TRUE only if vertical urban/land surface and trees exist)
+                surface_reflections = .TRUE.,         & !< flag to switch the calculation of radiation interaction between surfaces.
+                                                        !< When it switched off, only the effect of buildings and trees shadow
+                                                        !< will be considered. However fewer SVFs are expected.
+                radiation_interactions_on = .TRUE.      !< namelist flag to force RTM activiation regardless to vertical urban/land surface and trees
+    REAL(wp) :: albedo = 9999999.9_wp,           & !< NAMELIST alpha
+                albedo_lw_dif = 9999999.9_wp,    & !< NAMELIST aldif
+                albedo_lw_dir = 9999999.9_wp,    & !< NAMELIST aldir
+                albedo_sw_dif = 9999999.9_wp,    & !< NAMELIST asdif
+                albedo_sw_dir = 9999999.9_wp,    & !< NAMELIST asdir
+                decl_1,                          & !< declination coef. 1
+                decl_2,                          & !< declination coef. 2
+                decl_3,                          & !< declination coef. 3
+                dt_radiation = 0.0_wp,           & !< radiation model timestep
+                emissivity = 9999999.9_wp,       & !< NAMELIST surface emissivity
+                lon = 0.0_wp,                    & !< longitude in radians
+                lat = 0.0_wp,                    & !< latitude in radians
+                net_radiation = 0.0_wp,          & !< net radiation at surface
+                skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
+                sky_trans,                       & !< sky transmissivity
+                time_radiation = 0.0_wp,         & !< time since last call of radiation code
+                trace_fluxes_above = -1.0_wp,    & !< NAMELIST option for debug tracing of large radiative fluxes (W/m2;W/m3)
+                min_stable_coszen = 0.0262_wp      !< 1.5 deg above horizon, eliminates most of circumsolar
+    INTEGER(iwp) ::  day_of_year   !< day of the current year
+    REAL(wp) ::  cos_zenith        !< cosine of solar zenith angle, also z-coordinate of solar unit vector
+    REAL(wp) ::  d_hours_day       !< 1 / hours-per-day
+    REAL(wp) ::  d_seconds_hour    !< 1 / seconds-per-hour
+    REAL(wp) ::  second_of_day     !< second of the current day
+    REAL(wp) ::  sun_dir_lat       !< y-coordinate of solar unit vector
+    REAL(wp) ::  sun_dir_lon       !< x-coordinate of solar unit vector
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_net_av       !< average of net radiation (rad_net) at surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_in_xy_av  !< average of incoming longwave radiation at surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_out_xy_av !< average of outgoing longwave radiation at surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_in_xy_av  !< average of incoming shortwave radiation at surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_out_xy_av !< average of outgoing shortwave radiation at surface
+    REAL(wp), PARAMETER :: emissivity_atm_clsky = 0.8_wp       !< emissivity of the clear-sky atmosphere
+!-- Indices of radiation-related input attributes in building_surface_pars (others are in
+!-- urban_surface_mod)
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_b_wall  = 19  !< index for Broadband albedo of wall fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_l_wall  = 20  !< index for Longwave albedo of wall fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_s_wall  = 21  !< index for Shortwave albedo of wall fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_b_win   = 22  !< index for Broadband albedo of window fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_l_win   = 23  !< index for Longwave albedo of window fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_s_win   = 24  !< index for Shortwave albedo of window fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_b_green = 24  !< index for Broadband albedo of green fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_l_green = 25  !< index for Longwave albedo of green fraction
+    INTEGER(iwp), PARAMETER ::  ind_s_alb_s_green = 26  !< index for Shortwave albedo of green fraction
+    INTEGER(iwp) ::  albedo_type  = 9999999_iwp, &  !< Albedo surface type
+                     dots_rad     = 0_iwp           !< starting index for timeseries output
+    INTEGER(iwp) ::  day_of_year                    !< day of the current year
+    LOGICAL ::  unscheduled_radiation_calls = .TRUE., &  !< flag parameter indicating whether additional calls of the radiation
+                                                         !< code are allowed
+                constant_albedo =            .FALSE., &  !< flag parameter indicating whether the albedo may change depending on
+                                                         !< zenith
+                force_radiation_call =       .FALSE., &  !< flag parameter for unscheduled radiation calls
+                lw_radiation =                .TRUE., &  !< flag parameter indicating whether longwave radiation shall be calculated
+                radiation =                  .FALSE., &  !< flag parameter indicating whether the radiation model is used
+                sun_up =                      .TRUE., &  !< flag parameter indicating whether the sun is up or down
+                sw_radiation =                .TRUE., &  !< flag parameter indicating whether shortwave radiation shall be
+                                                         !< calculated
+                sun_direction =              .FALSE., &  !< flag parameter indicating whether solar direction shall be calculated
+                average_radiation =          .FALSE., &  !< flag to set the calculation of radiation averaging for the domain
+                radiation_interactions =     .FALSE., &  !< flag to activiate RTM (TRUE only if vertical urban/land surface and
+                                                         !< trees exist)
+                surface_reflections =         .TRUE., &  !< flag to switch the calculation of radiation interaction between
+                                                         !< surfaces. When it switched off, only the effect of buildings and trees
+                                                         !< shadow will be considered. However fewer SVFs are expected.
+                radiation_interactions_on = .TRUE.       !< namelist flag to force RTM activiation regardless to vertical urban/
+                                                         !< land surface and trees
+    REAL(wp) ::  albedo = 9999999.9_wp,           &  !< NAMELIST alpha
+                 albedo_lw_dif = 9999999.9_wp,    &  !< NAMELIST aldif
+                 albedo_lw_dir = 9999999.9_wp,    &  !< NAMELIST aldir
+                 albedo_sw_dif = 9999999.9_wp,    &  !< NAMELIST asdif
+                 albedo_sw_dir = 9999999.9_wp,    &  !< NAMELIST asdir
+                 decl_1,                          &  !< declination coef. 1
+                 decl_2,                          &  !< declination coef. 2
+                 decl_3,                          &  !< declination coef. 3
+                 dt_radiation = 0.0_wp,           &  !< radiation model timestep
+                 emissivity = 9999999.9_wp,       &  !< NAMELIST surface emissivity
+                 lon = 0.0_wp,                    &  !< longitude in radians
+                 lat = 0.0_wp,                    &  !< latitude in radians
+                 net_radiation = 0.0_wp,          &  !< net radiation at surface
+                 skip_time_do_radiation = 0.0_wp, &  !< Radiation model is not called before this time
+                 sky_trans,                       &  !< sky transmissivity
+                 time_radiation = 0.0_wp,         &  !< time since last call of radiation code
+                 trace_fluxes_above = -1.0_wp,    &  !< NAMELIST option for debug tracing of large radiative fluxes (W/m2;W/m3)
+                 min_stable_coszen = 0.0262_wp       !< 1.5 deg above horizon, eliminates most of circumsolar
+    REAL(wp) ::  cos_zenith      !< cosine of solar zenith angle, also z-coordinate of solar unit vector
+    REAL(wp) ::  d_hours_day     !< 1 / hours-per-day
+    REAL(wp) ::  d_seconds_hour  !< 1 / seconds-per-hour
+    REAL(wp) ::  second_of_day   !< second of the current day
+    REAL(wp) ::  sun_dir_lat     !< y-coordinate of solar unit vector
+    REAL(wp) ::  sun_dir_lon     !< x-coordinate of solar unit vector
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_net_av        !< average of net radiation (rad_net) at surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_in_xy_av   !< average of incoming longwave radiation at surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_out_xy_av  !< average of outgoing longwave radiation at surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_in_xy_av   !< average of incoming shortwave radiation at surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_out_xy_av  !< average of outgoing shortwave radiation at surface
+    REAL(wp), PARAMETER ::  emissivity_atm_clsky = 0.8_wp  !< emissivity of the clear-sky atmosphere
+!
 !-- Land surface albedos for solar zenith angle of 60degree after Briegleb (1992)
 !-- (broadband, longwave, shortwave ):   bb,      lw,      sw,
     REAL(wp), DIMENSION(0:2,1:33), PARAMETER :: albedo_pars = RESHAPE( (/&
 .06_wp, 0.06_wp, 0.06_wp,            & !  1
 .19_wp, 0.28_wp, 0.09_wp,            & !  2
 .23_wp, 0.33_wp, 0.11_wp,            & !  3
 .23_wp, 0.33_wp, 0.11_wp,            & !  4
 .25_wp, 0.34_wp, 0.14_wp,            & !  5
 .14_wp, 0.22_wp, 0.06_wp,            & !  6
 .17_wp, 0.27_wp, 0.06_wp,            & !  7
 .19_wp, 0.31_wp, 0.06_wp,            & !  8
 .14_wp, 0.22_wp, 0.06_wp,            & !  9
 .18_wp, 0.28_wp, 0.06_wp,            & ! 10
 .43_wp, 0.51_wp, 0.35_wp,            & ! 11
 .32_wp, 0.40_wp, 0.24_wp,            & ! 12
 .19_wp, 0.27_wp, 0.10_wp,            & ! 13
 .77_wp, 0.65_wp, 0.90_wp,            & ! 14
 .77_wp, 0.65_wp, 0.90_wp,            & ! 15
 .82_wp, 0.70_wp, 0.95_wp,            & ! 16
 .08_wp, 0.08_wp, 0.08_wp,            & ! 17
 .17_wp, 0.17_wp, 0.17_wp,            & ! 18
 .17_wp, 0.17_wp, 0.17_wp,            & ! 19
 .30_wp, 0.30_wp, 0.30_wp,            & ! 20
 .17_wp, 0.17_wp, 0.17_wp,            & ! 21
 .17_wp, 0.17_wp, 0.17_wp,            & ! 22
 .17_wp, 0.17_wp, 0.17_wp,            & ! 23
 .17_wp, 0.17_wp, 0.17_wp,            & ! 24
 .17_wp, 0.17_wp, 0.17_wp,            & ! 25
 .17_wp, 0.17_wp, 0.17_wp,            & ! 26
 .17_wp, 0.17_wp, 0.17_wp,            & ! 27
 .17_wp, 0.17_wp, 0.17_wp,            & ! 28
 .17_wp, 0.17_wp, 0.17_wp,            & ! 29
 .17_wp, 0.17_wp, 0.17_wp,            & ! 30
 .17_wp, 0.17_wp, 0.17_wp,            & ! 31
 .17_wp, 0.17_wp, 0.17_wp,            & ! 32
 .17_wp, 0.17_wp, 0.17_wp             & ! 33
                                  /), (/ 3, 33 /) )
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
                         rad_lw_cs_hr,                  & !< longwave clear sky radiation heating rate (K/s)
                         rad_lw_cs_hr_av,               & !< average of rad_lw_cs_hr
                         rad_lw_hr,                     & !< longwave radiation heating rate (K/s)
                         rad_lw_hr_av,                  & !< average of rad_sw_hr
                         rad_lw_in,                     & !< incoming longwave radiation (W/m2)
                         rad_lw_in_av,                  & !< average of rad_lw_in
                         rad_lw_out,                    & !< outgoing longwave radiation (W/m2)
                         rad_lw_out_av,                 & !< average of rad_lw_out
                         rad_sw_cs_hr,                  & !< shortwave clear sky radiation heating rate (K/s)
                         rad_sw_cs_hr_av,               & !< average of rad_sw_cs_hr
                         rad_sw_hr,                     & !< shortwave radiation heating rate (K/s)
                         rad_sw_hr_av,                  & !< average of rad_sw_hr
                         rad_sw_in,                     & !< incoming shortwave radiation (W/m2)
                         rad_sw_in_av,                  & !< average of rad_sw_in
                         rad_sw_out,                    & !< outgoing shortwave radiation (W/m2)
                         rad_sw_out_av                    !< average of rad_sw_out
+!-- (broadband, longwave, shortwave ):              bb,      lw,      sw,
+    REAL(wp), DIMENSION(0:2,1:33), PARAMETER ::  albedo_pars = RESHAPE( (/&
+.06_wp, 0.06_wp, 0.06_wp,                  &  !  1
+.19_wp, 0.28_wp, 0.09_wp,                  &  !  2
+.23_wp, 0.33_wp, 0.11_wp,                  &  !  3
+.23_wp, 0.33_wp, 0.11_wp,                  &  !  4
+.25_wp, 0.34_wp, 0.14_wp,                  &  !  5
+.14_wp, 0.22_wp, 0.06_wp,                  &  !  6
+.17_wp, 0.27_wp, 0.06_wp,                  &  !  7
+.19_wp, 0.31_wp, 0.06_wp,                  &  !  8
+.14_wp, 0.22_wp, 0.06_wp,                  &  !  9
+.18_wp, 0.28_wp, 0.06_wp,                  &  ! 10
+.43_wp, 0.51_wp, 0.35_wp,                  &  ! 11
+.32_wp, 0.40_wp, 0.24_wp,                  &  ! 12
+.19_wp, 0.27_wp, 0.10_wp,                  &  ! 13
+.77_wp, 0.65_wp, 0.90_wp,                  &  ! 14
+.77_wp, 0.65_wp, 0.90_wp,                  &  ! 15
+.82_wp, 0.70_wp, 0.95_wp,                  &  ! 16
+.08_wp, 0.08_wp, 0.08_wp,                  &  ! 17
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 18
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 19
+.30_wp, 0.30_wp, 0.30_wp,                  &  ! 20
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 21
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 22
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 23
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 24
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 25
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 26
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 27
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 28
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 29
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 30
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 31
+.17_wp, 0.17_wp, 0.17_wp,                  &  ! 32
+.17_wp, 0.17_wp, 0.17_wp                   &  ! 33
+                                               /), (/ 3, 33 /) )
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  &
+                        rad_lw_cs_hr,                   &  !< longwave clear sky radiation heating rate (K/s)
+                        rad_lw_cs_hr_av,                &  !< average of rad_lw_cs_hr
+                        rad_lw_hr,                      &  !< longwave radiation heating rate (K/s)
+                        rad_lw_hr_av,                   &  !< average of rad_sw_hr
+                        rad_lw_in,                      &  !< incoming longwave radiation (W/m2)
+                        rad_lw_in_av,                   &  !< average of rad_lw_in
+                        rad_lw_out,                     &  !< outgoing longwave radiation (W/m2)
+                        rad_lw_out_av,                  &  !< average of rad_lw_out
+                        rad_sw_cs_hr,                   &  !< shortwave clear sky radiation heating rate (K/s)
+                        rad_sw_cs_hr_av,                &  !< average of rad_sw_cs_hr
+                        rad_sw_hr,                      &  !< shortwave radiation heating rate (K/s)
+                        rad_sw_hr_av,                   &  !< average of rad_sw_hr
+                        rad_sw_in,                      &  !< incoming shortwave radiation (W/m2)
+                        rad_sw_in_av,                   &  !< average of rad_sw_in
+                        rad_sw_out,                     &  !< outgoing shortwave radiation (W/m2)
+                        rad_sw_out_av                      !< average of rad_sw_out
 …
 !-- Variables and parameters used in RRTMG only
 #if defined ( __rrtmg )
     CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
+    CHARACTER(LEN=12) ::  rrtm_input_file = "RAD_SND_DATA"  !< name of the NetCDF input file (sounding data)
+!
 !-- Flag parameters to be passed to RRTMG (should not be changed until ice phase in clouds is allowed)
     INTEGER(iwp), PARAMETER :: rrtm_idrv     = 1, & !< flag for longwave upward flux calculation option (0,1)
                                rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
                                rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
                                rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
                                rrtm_inflgsw  = 2, & !< flag for sw cloud optical properties (0,1,2)
                                rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
                                rrtm_liqflgsw = 1    !< flag for sw liquid droplet specifications
+!
 !-- The following variables should be only changed with care, as this will
 !-- require further setting of some variables, which is currently not
+!-- implemented (aerosols, ice phase).
     INTEGER(iwp) :: nzt_rad,           & !< upper vertical limit for radiation calculations
                     rrtm_icld = 0,     & !< cloud flag (0: clear sky column, 1: cloudy column)
                     rrtm_iaer = 0        !< aerosol option flag (0: no aerosol layers, for lw only: 6 (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented)
+    INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
     LOGICAL :: snd_exists = .FALSE.      !< flag parameter to check whether a user-defined input files exists
     LOGICAL :: sw_exists = .FALSE.       !< flag parameter to check whether that required rrtmg sw file exists
+    LOGICAL :: lw_exists = .FALSE.       !< flag parameter to check whether that required rrtmg lw file exists
+    REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
     REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !< hypostatic pressure from sounding data (hPa)
                                            rrtm_tsfc,   & !< dummy array for storing surface temperature
                                            t_snd          !< actual temperature from sounding data (hPa)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rrtm_ccl4vmr,   & !< CCL4 volume mixing ratio (g/mol)
                                              rrtm_cfc11vmr,  & !< CFC11 volume mixing ratio (g/mol)
                                              rrtm_cfc12vmr,  & !< CFC12 volume mixing ratio (g/mol)
                                              rrtm_cfc22vmr,  & !< CFC22 volume mixing ratio (g/mol)
                                              rrtm_ch4vmr,    & !< CH4 volume mixing ratio
                                              rrtm_cicewp,    & !< in-cloud ice water path (g/m2)
                                              rrtm_cldfr,     & !< cloud fraction (0,1)
                                              rrtm_cliqwp,    & !< in-cloud liquid water path (g/m2)
                                              rrtm_co2vmr,    & !< CO2 volume mixing ratio (g/mol)
                                              rrtm_emis,      & !< surface emissivity (0-1)
                                              rrtm_h2ovmr,    & !< H2O volume mixing ratio
                                              rrtm_n2ovmr,    & !< N2O volume mixing ratio
                                              rrtm_o2vmr,     & !< O2 volume mixing ratio
                                              rrtm_o3vmr,     & !< O3 volume mixing ratio
                                              rrtm_play,      & !< pressure layers (hPa, zu-grid)
                                              rrtm_plev,      & !< pressure layers (hPa, zw-grid)
                                              rrtm_reice,     & !< cloud ice effective radius (microns)
                                              rrtm_reliq,     & !< cloud water drop effective radius (microns)
                                              rrtm_tlay,      & !< actual temperature (K, zu-grid)
                                              rrtm_tlev,      & !< actual temperature (K, zw-grid)
                                              rrtm_lwdflx,    & !< RRTM output of incoming longwave radiation flux (W/m2)
                                              rrtm_lwdflxc,   & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
                                              rrtm_lwuflx,    & !< RRTM output of outgoing longwave radiation flux (W/m2)
                                              rrtm_lwuflxc,   & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
                                              rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
                                              rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
                                              rrtm_lwhr,      & !< RRTM output of longwave radiation heating rate (K/d)
                                              rrtm_lwhrc,     & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
                                              rrtm_swdflx,    & !< RRTM output of incoming shortwave radiation flux (W/m2)
                                              rrtm_swdflxc,   & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
                                              rrtm_swuflx,    & !< RRTM output of outgoing shortwave radiation flux (W/m2)
                                              rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
                                              rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
                                              rrtm_swhrc,     & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
                                              rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m2)
                                              rrtm_difdflux     !< RRTM output of incoming diffuse shortwave (W/m2)
     REAL(wp), DIMENSION(1) ::                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
                                              rrtm_aldir,     & !< surface albedo for longwave direct radiation
                                              rrtm_asdif,     & !< surface albedo for shortwave diffuse radiation
                                              rrtm_asdir        !< surface albedo for shortwave direct radiation
+    INTEGER(iwp), PARAMETER ::  rrtm_idrv     = 1, &  !< flag for longwave upward flux calculation option (0,1)
+                                rrtm_inflglw  = 2, &  !< flag for lw cloud optical properties (0,1,2)
+                                rrtm_iceflglw = 0, &  !< flag for lw ice particle specifications (0,1,2,3)
+                                rrtm_liqflglw = 1, &  !< flag for lw liquid droplet specifications
+                                rrtm_inflgsw  = 2, &  !< flag for sw cloud optical properties (0,1,2)
+                                rrtm_iceflgsw = 0, &  !< flag for sw ice particle specifications (0,1,2,3)
+                                rrtm_liqflgsw = 1     !< flag for sw liquid droplet specifications
+!
+!-- The following variables should be only changed with care, as this will require further setting
+!-- of some variables, which is currently not implemented (aerosols, ice phase).
+    INTEGER(iwp) :: nzt_rad,       &  !< upper vertical limit for radiation calculations
+                    rrtm_icld = 0, &  !< cloud flag (0: clear sky column, 1: cloudy column)
+                    rrtm_iaer = 0     !< aerosol option flag (0: no aerosol layers, for lw only: 6
+                                      !< (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented)
+    INTEGER(iwp) ::  nc_stat          !< local variable for storin the result of netCDF calls for error message handling
+    LOGICAL ::  snd_exists = .FALSE.  !< flag parameter to check whether a user-defined input files exists
+    LOGICAL ::  sw_exists = .FALSE.   !< flag parameter to check whether that required rrtmg sw file exists
+    LOGICAL ::  lw_exists = .FALSE.   !< flag parameter to check whether that required rrtmg lw file exists
+    REAL(wp), PARAMETER ::  mol_mass_air_d_wv = 1.607793_wp  !< molecular weight dry air / water vapor
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  hyp_snd,   &       !< hypostatic pressure from sounding data (hPa)
+                                              rrtm_tsfc, &       !< dummy array for storing surface temperature
+                                              t_snd              !< actual temperature from sounding data (hPa)
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rrtm_ccl4vmr,   &  !< CCL4 volume mixing ratio (g/mol)
+                                              rrtm_cfc11vmr,  &  !< CFC11 volume mixing ratio (g/mol)
+                                              rrtm_cfc12vmr,  &  !< CFC12 volume mixing ratio (g/mol)
+                                              rrtm_cfc22vmr,  &  !< CFC22 volume mixing ratio (g/mol)
+                                              rrtm_ch4vmr,    &  !< CH4 volume mixing ratio
+                                              rrtm_cicewp,    &  !< in-cloud ice water path (g/m2)
+                                              rrtm_cldfr,     &  !< cloud fraction (0,1)
+                                              rrtm_cliqwp,    &  !< in-cloud liquid water path (g/m2)
+                                              rrtm_co2vmr,    &  !< CO2 volume mixing ratio (g/mol)
+                                              rrtm_emis,      &  !< surface emissivity (0-1)
+                                              rrtm_h2ovmr,    &  !< H2O volume mixing ratio
+                                              rrtm_n2ovmr,    &  !< N2O volume mixing ratio
+                                              rrtm_o2vmr,     &  !< O2 volume mixing ratio
+                                              rrtm_o3vmr,     &  !< O3 volume mixing ratio
+                                              rrtm_play,      &  !< pressure layers (hPa, zu-grid)
+                                              rrtm_plev,      &  !< pressure layers (hPa, zw-grid)
+                                              rrtm_reice,     &  !< cloud ice effective radius (microns)
+                                              rrtm_reliq,     &  !< cloud water drop effective radius (microns)
+                                              rrtm_tlay,      &  !< actual temperature (K, zu-grid)
+                                              rrtm_tlev,      &  !< actual temperature (K, zw-grid)
+                                              rrtm_lwdflx,    &  !< RRTM output of incoming longwave radiation flux (W/m2)
+                                              rrtm_lwdflxc,   &  !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
+                                              rrtm_lwuflx,    &  !< RRTM output of outgoing longwave radiation flux (W/m2)
+                                              rrtm_lwuflxc,   &  !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
+                                              rrtm_lwuflx_dt, &  !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
+                                              rrtm_lwuflxc_dt,&  !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
+                                              rrtm_lwhr,      &  !< RRTM output of longwave radiation heating rate (K/d)
+                                              rrtm_lwhrc,     &  !< RRTM output of incoming longwave clear sky radiation heating
+                                                                 !< rate (K/d)
+                                              rrtm_swdflx,    &  !< RRTM output of incoming shortwave radiation flux (W/m2)
+                                              rrtm_swdflxc,   &  !< RRTM output of outgoing clear sky shortwave radiation flux(W/m2)
+                                              rrtm_swuflx,    &  !< RRTM output of outgoing shortwave radiation flux (W/m2)
+                                              rrtm_swuflxc,   &  !< RRTM output of incoming clear sky shortwave radiation flux(W/m2)
+                                              rrtm_swhr,      &  !< RRTM output of shortwave radiation heating rate (K/d)
+                                              rrtm_swhrc,     &  !< RRTM output of incoming shortwave clear sky radiation heating
+                                                                 !< rate (K/d)
+                                              rrtm_dirdflux,  &  !< RRTM output of incoming direct shortwave (W/m2)
+                                              rrtm_difdflux      !< RRTM output of incoming diffuse shortwave (W/m2)
+    REAL(wp), DIMENSION(1) ::  rrtm_aldif,     &  !< surface albedo for longwave diffuse radiation
+                               rrtm_aldir,     &  !< surface albedo for longwave direct radiation
+                               rrtm_asdif,     &  !< surface albedo for shortwave diffuse radiation
+                               rrtm_asdir         !< surface albedo for shortwave direct radiation
+!
 !-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   & !< incoming clear sky longwave radiation (W/m2) (not used)
                                                 rad_lw_cs_out,  & !< outgoing clear sky longwave radiation (W/m2) (not used)
                                                 rad_sw_cs_in,   & !< incoming clear sky shortwave radiation (W/m2) (not used)
                                                 rad_sw_cs_out,  & !< outgoing clear sky shortwave radiation (W/m2) (not used)
                                                 rrtm_lw_tauaer, & !< lw aerosol optical depth
                                                 rrtm_lw_taucld, & !< lw in-cloud optical depth
                                                 rrtm_sw_taucld, & !< sw in-cloud optical depth
                                                 rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
                                                 rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
                                                 rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
                                                 rrtm_sw_tauaer, & !< sw aerosol optical depth
                                                 rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
                                                 rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
                                                 rrtm_sw_ecaer     !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   &  !< incoming clear sky longwave radiation (W/m2) (not used)
+                                                rad_lw_cs_out,  &  !< outgoing clear sky longwave radiation (W/m2) (not used)
+                                                rad_sw_cs_in,   &  !< incoming clear sky shortwave radiation (W/m2) (not used)
+                                                rad_sw_cs_out,  &  !< outgoing clear sky shortwave radiation (W/m2) (not used)
+                                                rrtm_lw_tauaer, &  !< lw aerosol optical depth
+                                                rrtm_lw_taucld, &  !< lw in-cloud optical depth
+                                                rrtm_sw_taucld, &  !< sw in-cloud optical depth
+                                                rrtm_sw_ssacld, &  !< sw in-cloud single scattering albedo
+                                                rrtm_sw_asmcld, &  !< sw in-cloud asymmetry parameter
+                                                rrtm_sw_fsfcld, &  !< sw in-cloud forward scattering fraction
+                                                rrtm_sw_tauaer, &  !< sw aerosol optical depth
+                                                rrtm_sw_ssaaer, &  !< sw aerosol single scattering albedo
+                                                rrtm_sw_asmaer, &  !< sw aerosol asymmetry parameter
+                                                rrtm_sw_ecaer      !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
 #endif
+!
 !-- Parameters of urban and land surface models
+    INTEGER(iwp)                                   ::  nz_urban                           !< number of layers of urban surface (will be calculated)
+    INTEGER(iwp)                                   ::  nz_plant                           !< number of layers of plant canopy (will be calculated)
+    INTEGER(iwp)                                   ::  nz_urban_b                         !< bottom layer of urban surface (will be calculated)
+    INTEGER(iwp)                                   ::  nz_urban_t                         !< top layer of urban surface (will be calculated)
+    INTEGER(iwp)                                   ::  nz_plant_t                         !< top layer of plant canopy (will be calculated)
+!-- parameters of urban and land surface models
+    INTEGER(iwp), PARAMETER                        ::  nzut_free = 3                      !< number of free layers above top of of topography
+    INTEGER(iwp), PARAMETER                        ::  ndsvf = 2                          !< number of dimensions of real values in SVF
+    INTEGER(iwp), PARAMETER                        ::  idsvf = 2                          !< number of dimensions of integer values in SVF
+    INTEGER(iwp), PARAMETER                        ::  ndcsf = 1                          !< number of dimensions of real values in CSF
+    INTEGER(iwp), PARAMETER                        ::  idcsf = 2                          !< number of dimensions of integer values in CSF
+    INTEGER(iwp), PARAMETER                        ::  kdcsf = 4                          !< number of dimensions of integer values in CSF calculation array
+    INTEGER(iwp), PARAMETER                        ::  id = 1                             !< position of d-index in surfl and surf
+    INTEGER(iwp), PARAMETER                        ::  iz = 2                             !< position of k-index in surfl and surf
+    INTEGER(iwp), PARAMETER                        ::  iy = 3                             !< position of j-index in surfl and surf
+    INTEGER(iwp), PARAMETER                        ::  ix = 4                             !< position of i-index in surfl and surf
+    INTEGER(iwp), PARAMETER                        ::  im = 5                             !< position of surface m-index in surfl and surf
+    INTEGER(iwp), PARAMETER                        ::  nidx_surf = 5                      !< number of indices in surfl and surf
+    INTEGER(iwp), PARAMETER                        ::  nsurf_type = 10                    !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
+    INTEGER(iwp), PARAMETER                        ::  iup_u    = 0                       !< 0 - index of urban upward surface (ground or roof)
+    INTEGER(iwp), PARAMETER                        ::  idown_u  = 1                       !< 1 - index of urban downward surface (overhanging)
+    INTEGER(iwp), PARAMETER                        ::  inorth_u = 2                       !< 2 - index of urban northward facing wall
+    INTEGER(iwp), PARAMETER                        ::  isouth_u = 3                       !< 3 - index of urban southward facing wall
+    INTEGER(iwp), PARAMETER                        ::  ieast_u  = 4                       !< 4 - index of urban eastward facing wall
+    INTEGER(iwp), PARAMETER                        ::  iwest_u  = 5                       !< 5 - index of urban westward facing wall
+    INTEGER(iwp), PARAMETER                        ::  iup_l    = 6                       !< 6 - index of land upward surface (ground or roof)
+    INTEGER(iwp), PARAMETER                        ::  inorth_l = 7                       !< 7 - index of land northward facing wall
+    INTEGER(iwp), PARAMETER                        ::  isouth_l = 8                       !< 8 - index of land southward facing wall
+    INTEGER(iwp), PARAMETER                        ::  ieast_l  = 9                       !< 9 - index of land eastward facing wall
+    INTEGER(iwp), PARAMETER                        ::  iwest_l  = 10                      !< 10- index of land westward facing wall
+    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1/)   !< surface normal direction x indices
+    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/)   !< surface normal direction y indices
+    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/)   !< surface normal direction z indices
+    REAL(wp),     DIMENSION(0:nsurf_type)          :: facearea                            !< area of single face in respective
+                                                                                          !< direction (will be calc'd)
+!-- indices and sizes of urban and land surface models
+    INTEGER(iwp)                                   ::  startland        !< start index of block of land and roof surfaces
+    INTEGER(iwp)                                   ::  endland          !< end index of block of land and roof surfaces
+    INTEGER(iwp)                                   ::  nlands           !< number of land and roof surfaces in local processor
+    INTEGER(iwp)                                   ::  startwall        !< start index of block of wall surfaces
+    INTEGER(iwp)                                   ::  endwall          !< end index of block of wall surfaces
+    INTEGER(iwp)                                   ::  nwalls           !< number of wall surfaces in local processor
+!-- indices needed for RTM netcdf output subroutines
+    INTEGER(iwp), PARAMETER                        :: nd = 5
+    CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
+    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_u = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
+    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_l = (/ iup_l, isouth_l, inorth_l, iwest_l, ieast_l /)
+    INTEGER(iwp), DIMENSION(0:nd-1)                :: dirstart
+    INTEGER(iwp), DIMENSION(0:nd-1)                :: dirend
+!-- indices and sizes of urban and land surface models
+    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surfl            !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x, m]
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_linear     !< dtto (linearly allocated array)
+    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surf             !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x, m]
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surf_linear      !< dtto (linearly allocated array)
+    INTEGER(iwp)                                   ::  nsurfl           !< number of all surfaces in local processor
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  nsurfs           !< array of number of all surfaces in individual processors
+    INTEGER(iwp)                                   ::  nsurf            !< global number of surfaces in index array of surfaces (nsurf = proc nsurfs)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfstart        !< starts of blocks of surfaces for individual processors in array surf (indexed from 1)
+                                                                        !< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1]
+!-- block variables needed for calculation of the plant canopy model inside the urban surface model
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pct              !< top layer of the plant canopy
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pch              !< heights of the plant canopy
+    INTEGER(iwp)                                   ::  npcbl = 0        !< number of the plant canopy gridboxes in local processor
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pcbl             !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i]
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw          !< array of absorbed sw radiation for local plant canopy box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir       !< array of absorbed direct sw radiation for local plant canopy box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif       !< array of absorbed diffusion sw radiation for local plant canopy box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw          !< array of absorbed lw radiation for local plant canopy box
+!-- configuration parameters (they can be setup in PALM config)
+    LOGICAL                                        ::  raytrace_mpi_rma = .TRUE.          !< use MPI RMA to access LAD and gridsurf from remote processes during raytracing
+    LOGICAL                                        ::  rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for
+                                                                                          !< reflected radiation (as opposed to all mutually visible pairs)
+    LOGICAL                                        ::  plant_lw_interact = .TRUE.         !< whether plant canopy interacts with LW radiation (in addition to SW)
+    INTEGER(iwp)                                   ::  mrt_nlevels = 0                    !< number of vertical boxes above surface for which to calculate MRT
+    LOGICAL                                        ::  mrt_skip_roof = .TRUE.             !< do not calculate MRT above roof surfaces
+    LOGICAL                                        ::  mrt_include_sw = .TRUE.            !< should MRT calculation include SW radiation as well?
+    INTEGER(wp)                                    ::  mrt_geom = 1                       !< method for MRT direction weights simulating a sphere or a human body
+    REAL(wp), DIMENSION(2)                         ::  mrt_geom_params = (/ .12_wp, .88_wp /)   !< parameters for the selected method
+    INTEGER(iwp)                                   ::  nrefsteps = 3                      !< number of reflection steps to perform
+    REAL(wp), PARAMETER                            ::  ext_coef = 0.6_wp                  !< extinction coefficient (a.k.a. alpha)
+    INTEGER(iwp), PARAMETER                        ::  rad_version_len = 10               !< length of identification string of rad version
+    CHARACTER(rad_version_len), PARAMETER          ::  rad_version = 'RAD v. 3.0'         !< identification of version of binary svf and restart files
+    INTEGER(iwp)                                   ::  raytrace_discrete_elevs = 40       !< number of discretization steps for elevation (nadir to zenith)
+    INTEGER(iwp)                                   ::  raytrace_discrete_azims = 80       !< number of discretization steps for azimuth (out of 360 degrees)
+    REAL(wp)                                       ::  max_raytracing_dist = -999.0_wp    !< maximum distance for raytracing (in metres)
+    REAL(wp)                                       ::  min_irrf_value = 1e-6_wp           !< minimum potential irradiance factor value for raytracing
+    REAL(wp), DIMENSION(1:30)                      ::  svfnorm_report_thresh = 1e21_wp    !< thresholds of SVF normalization values to report
+    INTEGER(iwp)                                   ::  svfnorm_report_num                 !< number of SVF normalization thresholds to report
+!-- radiation related arrays to be used in radiation_interaction routine
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rad_sw_in_dir    !< direct sw radiation
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rad_sw_in_diff   !< diffusion sw radiation
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rad_lw_in_diff   !< diffusion lw radiation
+!-- parameters required for RRTMG lower boundary condition
+    REAL(wp)                   :: albedo_urb      !< albedo value retuned to RRTMG boundary cond.
+    REAL(wp)                   :: emissivity_urb  !< emissivity value retuned to RRTMG boundary cond.
+    REAL(wp)                   :: t_rad_urb       !< temperature value retuned to RRTMG boundary cond.
+!-- type for calculation of svf
+    INTEGER(iwp) ::  nz_urban    !< number of layers of urban surface (will be calculated)
+    INTEGER(iwp) ::  nz_plant    !< number of layers of plant canopy (will be calculated)
+    INTEGER(iwp) ::  nz_urban_b  !< bottom layer of urban surface (will be calculated)
+    INTEGER(iwp) ::  nz_urban_t  !< top layer of urban surface (will be calculated)
+    INTEGER(iwp) ::  nz_plant_t  !< top layer of plant canopy (will be calculated)
+!-- Parameters of urban and land surface models
+    INTEGER(iwp), PARAMETER ::  nzut_free = 3    !< number of free layers above top of of topography
+    INTEGER(iwp), PARAMETER ::  ndsvf = 2        !< number of dimensions of real values in SVF
+    INTEGER(iwp), PARAMETER ::  idsvf = 2        !< number of dimensions of integer values in SVF
+    INTEGER(iwp), PARAMETER ::  ndcsf = 1        !< number of dimensions of real values in CSF
+    INTEGER(iwp), PARAMETER ::  idcsf = 2        !< number of dimensions of integer values in CSF
+    INTEGER(iwp), PARAMETER ::  kdcsf = 4        !< number of dimensions of integer values in CSF calculation array
+    INTEGER(iwp), PARAMETER ::  id = 1           !< position of d-index in surfl and surf
+    INTEGER(iwp), PARAMETER ::  iz = 2           !< position of k-index in surfl and surf
+    INTEGER(iwp), PARAMETER ::  iy = 3           !< position of j-index in surfl and surf
+    INTEGER(iwp), PARAMETER ::  ix = 4           !< position of i-index in surfl and surf
+    INTEGER(iwp), PARAMETER ::  im = 5           !< position of surface m-index in surfl and surf
+    INTEGER(iwp), PARAMETER ::  nidx_surf = 5    !< number of indices in surfl and surf
+    INTEGER(iwp), PARAMETER ::  nsurf_type = 10  !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
+    INTEGER(iwp), PARAMETER ::  iup_u    = 0     !< 0 - index of urban upward surface (ground or roof)
+    INTEGER(iwp), PARAMETER ::  idown_u  = 1     !< 1 - index of urban downward surface (overhanging)
+    INTEGER(iwp), PARAMETER ::  inorth_u = 2     !< 2 - index of urban northward facing wall
+    INTEGER(iwp), PARAMETER ::  isouth_u = 3     !< 3 - index of urban southward facing wall
+    INTEGER(iwp), PARAMETER ::  ieast_u  = 4     !< 4 - index of urban eastward facing wall
+    INTEGER(iwp), PARAMETER ::  iwest_u  = 5     !< 5 - index of urban westward facing wall
+    INTEGER(iwp), PARAMETER ::  iup_l    = 6     !< 6 - index of land upward surface (ground or roof)
+    INTEGER(iwp), PARAMETER ::  inorth_l = 7     !< 7 - index of land northward facing wall
+    INTEGER(iwp), PARAMETER ::  isouth_l = 8     !< 8 - index of land southward facing wall
+    INTEGER(iwp), PARAMETER ::  ieast_l  = 9     !< 9 - index of land eastward facing wall
+    INTEGER(iwp), PARAMETER ::  iwest_l  = 10    !< 10- index of land westward facing wall
+    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1/)  !< surface normal direction x indic.
+    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/)  !< surface normal direction y indic.
+    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/)  !< surface normal direction z indic.
+    REAL(wp), DIMENSION(0:nsurf_type) ::  facearea  !< area of single face in respective direction (will be calc'd)
+!
+!-- Indices and sizes of urban and land surface models
+    INTEGER(iwp) ::  startland  !< start index of block of land and roof surfaces
+    INTEGER(iwp) ::  endland    !< end index of block of land and roof surfaces
+    INTEGER(iwp) ::  nlands     !< number of land and roof surfaces in local processor
+    INTEGER(iwp) ::  startwall  !< start index of block of wall surfaces
+    INTEGER(iwp) ::  endwall    !< end index of block of wall surfaces
+    INTEGER(iwp) ::  nwalls     !< number of wall surfaces in local processor
+!
+!-- Indices needed for RTM netcdf output subroutines
+    INTEGER(iwp), PARAMETER                        ::  nd = 5  !<
+    CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER ::  dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)  !<
+    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     ::  dirint_u = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)  !<
+    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     ::  dirint_l = (/ iup_l, isouth_l, inorth_l, iwest_l, ieast_l /)  !<
+    INTEGER(iwp), DIMENSION(0:nd-1)                ::  dirstart  !<
+    INTEGER(iwp), DIMENSION(0:nd-1)                ::  dirend    !<
+!
+!-- Indices and sizes of urban and land surface models
+    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surfl         !< coordinates of i-th local surface in local grid
+                                                                     !< - surfl[:,k] = [d, z, y, x, m]
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_linear  !< dtto (linearly allocated array)
+    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surf          !< coordinates of i-th surface in grid
+                                                                     !< - surf[:,k] = [d, z, y, x, m]
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surf_linear   !< dtto (linearly allocated array)
+    INTEGER(iwp)                                   ::  nsurfl        !< number of all surfaces in local processor
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  nsurfs        !< array of number of all surfaces in individual processors
+    INTEGER(iwp)                                   ::  nsurf         !< global number of surfaces in index array of surfaces
+                                                                     !< (nsurf = proc nsurfs)
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfstart     !< starts of blocks of surfaces for individual processors in
+                                                                     !< array surf (indexed from 1)
+                          !< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1]
+!
+!-- Block variables needed for calculation of the plant canopy model inside the urban surface model
+    INTEGER(iwp)                              ::  npcbl = 0   !< number of the plant canopy gridboxes in local processor
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pct         !< top layer of the plant canopy
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pch         !< heights of the plant canopy
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pcbl        !< k,j,i coordinates of l-th local plant canopy box
+                                                              !< pcbl[:,l] = [k, j, i]
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  pcbinsw     !< array of absorbed sw radiation for local plant canopy box
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  pcbinswdir  !< array of absorbed direct sw radiation for local plant canopy box
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  pcbinswdif  !< array of absorbed diffusion sw radiation for local plant canopy box
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  pcbinlw     !< array of absorbed lw radiation for local plant canopy box
+!
+!-- Configuration parameters (they can be setup in PALM config)
+    INTEGER(iwp), PARAMETER ::  rad_version_len = 10  !< length of identification string of rad version
+    CHARACTER(rad_version_len), PARAMETER ::  rad_version = 'RAD v. 3.0'  !< identification of version of binary svf and restart
+                                                                          !< files
+    INTEGER(iwp) ::  mrt_nlevels = 0               !< number of vertical boxes above surface for which to calculate MRT
+    INTEGER(iwp) ::  svfnorm_report_num            !< number of SVF normalization thresholds to report
+    INTEGER(iwp) ::  raytrace_discrete_elevs = 40  !< number of discretization steps for elevation (nadir to zenith)
+    INTEGER(iwp) ::  raytrace_discrete_azims = 80  !< number of discretization steps for azimuth (out of 360 degrees)
+    INTEGER(wp)  ::  mrt_geom = 1                  !< method for MRT direction weights simulating a sphere or a human body
+    INTEGER(iwp) ::  nrefsteps = 3                 !< number of reflection steps to perform
+    LOGICAL ::  raytrace_mpi_rma = .TRUE.            !< use MPI RMA to access LAD and gridsurf from remote processes during
+                                                     !< raytracing
+    LOGICAL ::  rad_angular_discretization = .TRUE.  !< whether to use fixed resolution discretization of view factors for
+                                                     !< reflected radiation (as opposed to all mutually visible pairs)
+    LOGICAL ::  plant_lw_interact = .TRUE.           !< whether plant canopy interacts with LW radiation (in addition to SW)
+    LOGICAL ::  mrt_skip_roof = .TRUE.               !< do not calculate MRT above roof surfaces
+    LOGICAL ::  mrt_include_sw = .TRUE.              !< should MRT calculation include SW radiation as well?
+    REAL(wp) ::  max_raytracing_dist = -999.0_wp     !< maximum distance for raytracing (in metres)
+    REAL(wp) ::  min_irrf_value = 1.0E-6_wp          !< minimum potential irradiance factor value for raytracing
+    REAL(wp), PARAMETER ::  ext_coef = 0.6_wp  !< extinction coefficient (a.k.a. alpha)
+    REAL(wp), DIMENSION(2)    ::  mrt_geom_params = (/ .12_wp, .88_wp /)  !< parameters for the selected method
+    REAL(wp), DIMENSION(1:30) ::  svfnorm_report_thresh = 1e21_wp         !< thresholds of SVF normalization values to report
+!
+!-- Radiation related arrays to be used in radiation_interaction routine
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_in_dir   !< direct sw radiation
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_in_diff  !< diffusion sw radiation
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_in_diff  !< diffusion lw radiation
+!
+!-- Parameters required for RRTMG lower boundary condition
+    REAL(wp) :: albedo_urb      !< albedo value retuned to RRTMG boundary cond.
+    REAL(wp) :: emissivity_urb  !< emissivity value retuned to RRTMG boundary cond.
+    REAL(wp) :: t_rad_urb       !< temperature value retuned to RRTMG boundary cond.
+!
+!-- Type for calculation of svf
     TYPE t_svf
         INTEGER(iwp)                               :: isurflt           !<
         INTEGER(iwp)                               :: isurfs            !<
         REAL(wp)                                   :: rsvf              !<
         REAL(wp)                                   :: rtransp           !<
+        INTEGER(iwp) ::  isurflt  !<
+        INTEGER(iwp) ::  isurfs   !<
+        REAL(wp)     ::  rsvf     !<
+        REAL(wp)     ::  rtransp  !<
     END TYPE
 !-- type for calculation of csf
+!
+!-- Type for calculation of csf
     TYPE t_csf
+        INTEGER(iwp)                               :: ip                !<
+        INTEGER(iwp)                               :: itx               !<
+        INTEGER(iwp)                               :: ity               !<
+        INTEGER(iwp)                               :: itz               !<
+        INTEGER(iwp)                               :: isurfs            !< Idx of source face / -1 for sky
+        REAL(wp)                                   :: rcvf              !< Canopy view factor for faces /
+                                                                        !< canopy sink factor for sky (-1)
+        INTEGER(iwp) ::  ip      !<
+        INTEGER(iwp) ::  itx     !<
+        INTEGER(iwp) ::  ity     !<
+        INTEGER(iwp) ::  itz     !<
+        INTEGER(iwp) ::  isurfs  !< Idx of source face / -1 for sky
+        REAL(wp)     ::  rcvf    !< Canopy view factor for faces / canopy sink factor for sky (-1)
     END TYPE
+!-- arrays storing the values of USM
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  svfsurf          !< svfsurf[:,isvf] = index of target and source surface for svf[isvf]
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  svf              !< array of shape view factors+direct irradiation factors for local surfaces
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins          !< array of sw radiation falling to local surface after i-th reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl          !< array of lw radiation for local surface after i-th reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  skyvf            !< array of sky view factor for each local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  skyvft           !< array of sky view factor including transparency for each local surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsitrans         !< dsidir[isvfl,i] = path transmittance of i-th
+                                                                        !< direction of direct solar irradiance per target surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsitransc        !< dtto per plant canopy box
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsidir           !< dsidir[:,i] = unit vector of i-th
+                                                                        !< direction of direct solar irradiance
+    INTEGER(iwp)                                   ::  ndsidir          !< number of apparent solar directions used
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  dsidir_rev       !< dsidir_rev[ielev,iazim] = i for dsidir or -1 if not present
+    INTEGER(iwp)                                   ::  nmrtbl           !< No. of local grid boxes for which MRT is calculated
+    INTEGER(iwp)                                   ::  nmrtf            !< number of MRT factors for local processor
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  mrtbl            !< coordinates of i-th local MRT box - surfl[:,i] = [z, y, x]
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  mrtfsurf         !< mrtfsurf[:,imrtf] = index of target MRT box and source surface for mrtf[imrtf]
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtf             !< array of MRT factors for each local MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtft            !< array of MRT factors including transparency for each local MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtsky           !< array of sky view factor for each local MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtskyt          !< array of sky view factor including transparency for each local MRT box
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  mrtdsit          !< array of direct solar transparencies for each local MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinsw          !< mean SW radiant flux for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinlw          !< mean LW radiant flux for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrt              !< mean radiant temperature for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinsw_av       !< time average mean SW radiant flux for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinlw_av       !< time average mean LW radiant flux for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrt_av           !< time average mean radiant temperature for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw         !< array of sw radiation falling to local surface including radiation from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw         !< array of lw radiation falling to local surface including radiation from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir      !< array of direct sw radiation falling to local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif      !< array of diffuse sw radiation from sky and model boundary falling to local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif      !< array of diffuse lw radiation from sky and model boundary falling to local surface
+                                                                        !< Outward radiation is only valid for nonvirtual surfaces
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsl        !< array of reflected sw radiation for local surface in i-th reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutll        !< array of reflected + emitted lw radiation for local surface in i-th reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfouts         !< array of reflected sw radiation for all surfaces in i-th reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutl         !< array of reflected + emitted lw radiation for all surfaces in i-th reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlg         !< global array of incoming lw radiation from plant canopy
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw        !< array of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw        !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfemitlwl      !< array of emitted lw radiation for local surface used to calculate effective surface temperature for radiation model
+!-- block variables needed for calculation of the plant canopy model inside the urban surface model
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  csfsurf          !< csfsurf[:,icsf] = index of target surface and csf grid index for csf[icsf]
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  csf              !< array of plant canopy sink fators + direct irradiation factors (transparency)
+    REAL(wp), DIMENSION(:,:,:), POINTER            ::  sub_lad          !< subset of lad_s within urban surface, transformed to plain Z coordinate
+!
+!-- Arrays storing the values of USM
+    INTEGER(iwp)                              ::  ndsidir      !< number of apparent solar directions used
+    INTEGER(iwp)                              ::  nmrtbl       !< No. of local grid boxes for which MRT is calculated
+    INTEGER(iwp)                              ::  nmrtf        !< number of MRT factors for local processor
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  mrtbl        !< coordinates of i-th local MRT box - surfl[:,i] = [z, y, x]
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  mrtfsurf     !< mrtfsurf[:,imrtf] = index of target MRT box and source surface for
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  dsidir_rev   !< dsidir_rev[ielev,iazim] = i for dsidir or -1 if not present
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  svfsurf      !< svfsurf[:,isvf] = index of target and source surface for svf[isvf]
+                                                               !< mrtf[imrtf]
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfins      !< array of sw radiation falling to local surface after i-th
+                                                               !< reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinl      !< array of lw radiation for local surface after i-th reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  skyvf        !< array of sky view factor for each local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  skyvft       !< array of sky view factor including transparency for each local
+                                                               !< surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtf         !< array of MRT factors for each local MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtft        !< array of MRT factors including transparency for each local MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtsky       !< array of sky view factor for each local MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtskyt      !< array of sky view factor including transparency for each local
+                                                               !< MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtinsw      !< mean SW radiant flux for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtinlw      !< mean LW radiant flux for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrt          !< mean radiant temperature for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtinsw_av   !< time average mean SW radiant flux for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtinlw_av   !< time average mean LW radiant flux for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrt_av       !< time average mean radiant temperature for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinsw     !< array of sw radiation falling to local surface including radiation
+                                                               !< from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinlw     !< array of lw radiation falling to local surface including radiation
+                                                               !< from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinswdir  !< array of direct sw radiation falling to local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinswdif  !< array of diffuse sw radiation from sky and model boundary falling
+                                                               !< to local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinlwdif  !< array of diffuse lw radiation from sky and model boundary falling
+                                                               !< to local surface
+                                                               !< Outward radiation is only valid for nonvirtual surfaces
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfoutsl    !< array of reflected sw radiation for local surface in i-th
+                                                               !< reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfoutll    !< array of reflected + emitted lw radiation for local surface in
+                                                               !< i-th reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfouts     !< array of reflected sw radiation for all surfaces in i-th
+                                                               !< reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfoutl     !< array of reflected + emitted lw radiation for all surfaces in
+                                                               !< i-th reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinlg     !< global array of incoming lw radiation from plant canopy
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfoutsw    !< array of total sw radiation outgoing from nonvirtual surfaces
+                                                               !< surfaces after all reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfoutlw    !< array of total lw radiation outgoing from nonvirtual surfaces
+                                                               !< surfaces after all reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfemitlwl  !< array of emitted lw radiation for local surface used to calculate
+                                                               !< effective surface temperature for radiation model
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  svf          !< array of shape view factors+direct irradiation factors for local
+                                                               !< surfaces
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  mrtdsit      !< array of direct solar transparencies for each local MRT box
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  dsitrans     !< dsidir[isvfl,i] = path transmittance of i-th
+                                                               !< direction of direct solar irradiance per target surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  dsitransc    !< dtto per plant canopy box
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  dsidir       !< dsidir[:,i] = unit vector of i-th
+                                                               !< direction of direct solar irradiance
+!
+!-- Block variables needed for calculation of the plant canopy model inside the urban surface model
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  csfsurf  !< csfsurf[:,icsf] = index of target surface and csf grid index
+                                                           !< for csf[icsf]
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  csf      !< array of plant canopy sink fators + direct irradiation factors
+                                                           !< (transparency)
+    REAL(wp), DIMENSION(:,:,:), POINTER       ::  sub_lad  !< subset of lad_s within urban surface, transformed to plain
+                                                           !< Z coordinate
 #if defined( __parallel )
     REAL(wp), DIMENSION(:), POINTER                ::  sub_lad_g        !< sub_lad globalized (used to avoid MPI RMA calls in raytracing)
+    REAL(wp), DIMENSION(:), POINTER ::  sub_lad_g  !< sub_lad globalized (used to avoid MPI RMA calls in raytracing)
 #endif
+    REAL(wp)                                       ::  prototype_lad    !< prototype leaf area density for computing effective optical depth
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nzterr, plantt   !< temporary global arrays for raytracing
+    INTEGER(iwp)                                   ::  plantt_max
+!-- arrays and variables for calculation of svf and csf
+    TYPE(t_svf), DIMENSION(:), POINTER             ::  asvf             !< pointer to growing svc array
+    TYPE(t_csf), DIMENSION(:), POINTER             ::  acsf             !< pointer to growing csf array
+    TYPE(t_svf), DIMENSION(:), POINTER             ::  amrtf            !< pointer to growing mrtf array
+    TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  asvf1, asvf2     !< realizations of svf array
+    TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET ::  acsf1, acsf2     !< realizations of csf array
+    TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  amrtf1, amrtf2   !< realizations of mftf array
+    INTEGER(iwp)                                   ::  nsvfla           !< dimmension of array allocated for storage of svf in local processor
+    INTEGER(iwp)                                   ::  ncsfla           !< dimmension of array allocated for storage of csf in local processor
+    INTEGER(iwp)                                   ::  nmrtfa           !< dimmension of array allocated for storage of mrt
+    INTEGER(iwp)                                   ::  msvf, mcsf, mmrtf!< mod for swapping the growing array
+    INTEGER(iwp), PARAMETER                        ::  gasize = 100000_iwp  !< initial size of growing arrays
+    REAL(wp), PARAMETER                            ::  grow_factor = 1.4_wp !< growth factor of growing arrays
+    INTEGER(iwp)                                   ::  nsvfl            !< number of svf for local processor
+    INTEGER(iwp)                                   ::  ncsfl            !< no. of csf in local processor
+                                                                        !< needed only during calc_svf but must be here because it is
+                                                                        !< shared between subroutines calc_svf and raytrace
+    INTEGER(iwp), DIMENSION(:,:,:,:), POINTER      ::  gridsurf         !< reverse index of local surfl[d,k,j,i] (for case rad_angular_discretization)
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE    ::  gridpcbl         !< reverse index of local pcbl[k,j,i]
+    INTEGER(iwp), PARAMETER                        ::  nsurf_type_u = 6 !< number of urban surf types (used in gridsurf)
+!-- temporary arrays for calculation of csf in raytracing
+    INTEGER(iwp)                                   ::  maxboxesg        !< max number of boxes ray can cross in the domain
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  boxes            !< coordinates of gridboxes being crossed by ray
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  crlens           !< array of crossing lengths of ray for particular grid boxes
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  lad_ip           !< array of numbers of process where lad is stored
+    INTEGER(iwp)                            ::  plantt_max
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nzterr, plantt  !< temporary global arrays for raytracing
+    REAL(wp)                                ::  prototype_lad   !< prototype leaf area density for computing effective optical depth
+!
+!-- Arrays and variables for calculation of svf and csf
+    TYPE(t_svf), DIMENSION(:), POINTER             ::  asvf                  !< pointer to growing svc array
+    TYPE(t_csf), DIMENSION(:), POINTER             ::  acsf                  !< pointer to growing csf array
+    TYPE(t_svf), DIMENSION(:), POINTER             ::  amrtf                 !< pointer to growing mrtf array
+    TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  asvf1, asvf2          !< realizations of svf array
+    TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET ::  acsf1, acsf2          !< realizations of csf array
+    TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  amrtf1, amrtf2        !< realizations of mftf array
+    INTEGER(iwp) ::  nsvfla             !< dimmension of array allocated for storage of svf in local processor
+    INTEGER(iwp) ::  ncsfla             !< dimmension of array allocated for storage of csf in local processor
+    INTEGER(iwp) ::  nmrtfa             !< dimmension of array allocated for storage of mrt
+    INTEGER(iwp) ::  msvf, mcsf, mmrtf  !< mod for swapping the growing array
+    INTEGER(iwp) ::  nsvfl              !< number of svf for local processor
+    INTEGER(iwp) ::  ncsfl              !< no. of csf in local processor needed only during calc_svf but must be here because
+                                        !< it is shared between subroutines calc_svf and raytrace
+    INTEGER(iwp), PARAMETER ::  gasize = 100000_iwp   !< initial size of growing arrays
+    INTEGER(iwp), PARAMETER ::  nsurf_type_u = 6      !< number of urban surf types (used in gridsurf)
+    INTEGER(iwp), DIMENSION(:,:,:,:), POINTER   ::  gridsurf              !< reverse index of local surfl[d,k,j,i] (for case
+                                                                          !< rad_angular_discretization)
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  gridpcbl              !< reverse index of local pcbl[k,j,i]
+    REAL(wp), PARAMETER ::  grow_factor = 1.4_wp  !< growth factor of growing arrays
+!
+!-- Temporary arrays for calculation of csf in raytracing
+    INTEGER(iwp) ::  maxboxesg  !< max number of boxes ray can cross in the domain
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  boxes      !< coordinates of gridboxes being crossed by ray
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  lad_ip     !< array of numbers of process where lad is stored
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  crlens     !< array of crossing lengths of ray for particular grid boxes
 #if defined( __parallel )
+    INTEGER(kind=MPI_ADDRESS_KIND), &
+                  DIMENSION(:), ALLOCATABLE        ::  lad_disp         !< array of displaycements of lad in local array of proc lad_ip
+    INTEGER(iwp)                                   ::  win_lad          !< MPI RMA window for leaf area density
+    INTEGER(iwp)                                   ::  win_gridsurf     !< MPI RMA window for reverse grid surface index
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  lad_s_ray        !< array of received lad_s for appropriate gridboxes crossed by ray
+    INTEGER(iwp) ::  win_lad       !< MPI RMA window for leaf area density
+    INTEGER(iwp) ::  win_gridsurf  !< MPI RMA window for reverse grid surface index
+    INTEGER(KIND=MPI_ADDRESS_KIND), DIMENSION(:), ALLOCATABLE ::  lad_disp  !< array of displaycements of lad in local array of
+                                                                            !< proc lad_ip
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad_s_ray  !< array of received lad_s for appropriate gridboxes crossed by ray
 #endif
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  target_surfl
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  rt2_track
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rt2_track_lad
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rt2_track_dist
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rt2_dist
+!-- arrays for time averages
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfradnet_av    !< average of net radiation to local surface including radiation from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw_av      !< average of sw radiation falling to local surface including radiation from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw_av      !< average of lw radiation falling to local surface including radiation from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir_av   !< average of direct sw radiation falling to local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif_av   !< average of diffuse sw radiation from sky and model boundary falling to local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif_av   !< average of diffuse lw radiation from sky and model boundary falling to local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswref_av   !< average of sw radiation falling to surface from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwref_av   !< average of lw radiation falling to surface from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw_av     !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw_av     !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw_av       !< Average of pcbinlw
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw_av       !< Average of pcbinsw
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir_av    !< Average of pcbinswdir
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif_av    !< Average of pcbinswdif
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswref_av    !< Average of pcbinswref
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  target_surfl    !<
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  rt2_track       !<
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rt2_track_dist  !<
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rt2_dist        !<
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rt2_track_lad   !<
+!
+!-- Arrays for time averages
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfradnet_av   !< average of net radiation to local surface including radiation from
+                                                            !< reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinsw_av     !< average of sw radiation falling to local surface including radiation
+                                                            !< from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlw_av     !< average of lw radiation falling to local surface including radiation
+                                                            !< from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinswdir_av  !< average of direct sw radiation falling to local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinswdif_av  !< average of diffuse sw radiation from sky and model boundary falling
+                                                            !< to local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlwdif_av  !< average of diffuse lw radiation from sky and model boundary falling
+                                                            !< to local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinswref_av  !< average of sw radiation falling to surface from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlwref_av  !< average of lw radiation falling to surface from reflections
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutsw_av    !< average of total sw radiation outgoing from nonvirtual surfaces
+                                                            !< surfaces after all reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutlw_av    !< average of total lw radiation outgoing from nonvirtual surfaces
+                                                            !< surfaces after all reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfins_av      !< average of array of residua of sw radiation absorbed in surface after
+                                                            !< last reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinl_av      !< average of array of residua of lw radiation absorbed in surface after
+                                                            !< last reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pcbinlw_av      !< Average of pcbinlw
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pcbinsw_av      !< Average of pcbinsw
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pcbinswdir_av   !< Average of pcbinswdir
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pcbinswdif_av   !< Average of pcbinswdif
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pcbinswref_av   !< Average of pcbinswref
 …
 !-- Energy balance variables
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !-- parameters of the land, roof and wall surfaces
     REAL(wp), DIMENSION(:), ALLOCATABLE            :: albedo_surf        !< albedo of the surface
     REAL(wp), DIMENSION(:), ALLOCATABLE            :: emiss_surf         !< emissivity of the wall surface
+!
 !-- External radiation. Depending on the given level of detail either a 1D or
 !-- a 3D array will be allocated.
     TYPE( real_1d_3d ) ::  rad_lw_in_f     !< external incoming longwave radiation, from observation or model
     TYPE( real_1d_3d ) ::  rad_sw_in_f     !< external incoming shortwave radiation, from observation or model
     TYPE( real_1d_3d ) ::  rad_sw_in_dif_f !< external incoming shortwave radiation, diffuse part, from observation or model
     TYPE( real_1d_3d ) ::  time_rad_f      !< time dimension for external radiation, from observation or model
+!-- Parameters of the land, roof and wall surfaces
+    REAL(wp), DIMENSION(:), ALLOCATABLE :: albedo_surf  !< albedo of the surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE :: emiss_surf   !< emissivity of the wall surface
+!
+!-- External radiation. Depending on the given level of detail either a 1D or a 3D array will be
+!-- allocated.
+    TYPE( real_1d_3d ) ::  rad_lw_in_f      !< external incoming longwave radiation, from observation or model
+    TYPE( real_1d_3d ) ::  rad_sw_in_f      !< external incoming shortwave radiation, from observation or model
+    TYPE( real_1d_3d ) ::  rad_sw_in_dif_f  !< external incoming shortwave radiation, diffuse part, from observation or model
+    TYPE( real_1d_3d ) ::  time_rad_f       !< time dimension for external radiation, from observation or model
     INTERFACE radiation_check_data_output
 …
+!
 !-- Public functions / NEEDS SORTING
+    PUBLIC radiation_check_data_output, radiation_check_data_output_pr,        &
+           radiation_check_data_output_ts,                                     &
+           radiation_check_parameters, radiation_control,                      &
+           radiation_header, radiation_init, radiation_parin,                  &
+           radiation_3d_data_averaging,                                        &
+           radiation_data_output_2d, radiation_data_output_3d,                 &
+           radiation_define_netcdf_grid, radiation_wrd_local,                  &
+           radiation_rrd_local, radiation_data_output_mask,                    &
+           radiation_calc_svf, radiation_write_svf,                            &
+           radiation_interaction, radiation_interaction_init,                  &
+           radiation_read_svf, radiation_presimulate_solar_pos
+    PUBLIC radiation_check_data_output,                                                            &
+           radiation_check_data_output_pr,                                                         &
+           radiation_check_data_output_ts,                                                         &
+           radiation_check_parameters,                                                             &
+           radiation_control,                                                                      &
+           radiation_header,                                                                       &
+           radiation_init,                                                                         &
+           radiation_parin,                                                                        &
+           radiation_3d_data_averaging,                                                            &
+           radiation_data_output_2d,                                                               &
+           radiation_data_output_3d,                                                               &
+           radiation_define_netcdf_grid,                                                           &
+           radiation_wrd_local,                                                                    &
+           radiation_rrd_local,                                                                    &
+           radiation_data_output_mask,                                                             &
+           radiation_calc_svf,                                                                     &
+           radiation_write_svf,                                                                    &
+           radiation_interaction,                                                                  &
+           radiation_interaction_init,                                                             &
+           radiation_read_svf,                                                                     &
+           radiation_presimulate_solar_pos
+!
 !-- Public variables and constants / NEEDS SORTING
+    PUBLIC albedo, albedo_type, decl_1, decl_2, decl_3, dots_rad, dt_radiation,&
+           emissivity, force_radiation_call, lat, lon, mrt_geom,               &
+           mrt_geom_params,                                                    &
+           mrt_include_sw, mrt_nlevels, mrtbl, mrtinsw, mrtinlw, nmrtbl,       &
+           rad_net_av, radiation, radiation_scheme, rad_lw_in,                 &
+           rad_lw_in_av, rad_lw_out, rad_lw_out_av,                            &
+           rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
+           rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
+           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, solar_constant,           &
+           skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
+           cos_zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,   &
+           idir, jdir, kdir, id, iz, iy, ix,                                   &
+           iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
+           iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
+           nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf,                 &
+           idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
+           radiation_interactions, startwall, startland, endland, endwall,     &
+           skyvf, skyvft, radiation_interactions_on, average_radiation,        &
+           rad_sw_in_diff, rad_sw_in_dir
+    PUBLIC albedo,                                                                                 &
+           albedo_type,                                                                            &
+           average_radiation,                                                                      &
+           calc_zenith,                                                                            &
+           cos_zenith,                                                                             &
+           decl_1,                                                                                 &
+           decl_2,                                                                                 &
+           decl_3,                                                                                 &
+           dots_rad,                                                                               &
+           dt_radiation,                                                                           &
+           endland,                                                                                &
+           endwall,                                                                                &
+           emissivity,                                                                             &
+           force_radiation_call,                                                                   &
+           id,                                                                                     &
+           idcsf,                                                                                  &
+           idir,                                                                                   &
+           idsvf,                                                                                  &
+           ieast_l,                                                                                &
+           ieast_u,                                                                                &
+           inorth_l,                                                                               &
+           inorth_u,                                                                               &
+           isouth_l,                                                                               &
+           isouth_u,                                                                               &
+           iup_l,                                                                                  &
+           iup_u,                                                                                  &
+           iwest_l,                                                                                &
+           iwest_u,                                                                                &
+           ix,                                                                                     &
+           iy,                                                                                     &
+           iz,                                                                                     &
+           jdir,                                                                                   &
+           kdcsf,                                                                                  &
+           kdir,                                                                                   &
+           lat,                                                                                    &
+           lon,                                                                                    &
+           mrtbl,                                                                                  &
+           mrt_geom,                                                                               &
+           mrt_geom_params,                                                                        &
+           mrt_include_sw,                                                                         &
+           mrt_nlevels,                                                                            &
+           mrtinsw,                                                                                &
+           mrtinlw,                                                                                &
+           nmrtbl,                                                                                 &
+           nsurf_type,                                                                             &
+           nz_urban_b,                                                                             &
+           nz_urban_t,                                                                             &
+           nz_urban,                                                                               &
+           nsurf,                                                                                  &
+           ndsvf,                                                                                  &
+           ndcsf,                                                                                  &
+           pch,                                                                                    &
+           pct,                                                                                    &
+           rad_net_av,                                                                             &
+           radiation,                                                                              &
+           radiation_scheme,                                                                       &
+           rad_lw_in,                                                                              &
+           rad_lw_in_av,                                                                           &
+           rad_lw_out,                                                                             &
+           rad_lw_out_av,                                                                          &
+           rad_lw_cs_hr,                                                                           &
+           rad_lw_cs_hr_av,                                                                        &
+           rad_lw_hr,                                                                              &
+           rad_lw_hr_av,                                                                           &
+           rad_sw_in,                                                                              &
+           rad_sw_in_av,                                                                           &
+           rad_sw_out,                                                                             &
+           rad_sw_out_av,                                                                          &
+           rad_sw_cs_hr,                                                                           &
+           rad_sw_cs_hr_av,                                                                        &
+           rad_sw_hr,                                                                              &
+           rad_sw_hr_av,                                                                           &
+           radiation_interactions,                                                                 &
+           radiation_interactions_on,                                                              &
+           rad_sw_in_diff,                                                                         &
+           rad_sw_in_dir,                                                                          &
+           solar_constant,                                                                         &
+           skip_time_do_radiation,                                                                 &
+           sun_direction,                                                                          &
+           sun_dir_lat,                                                                            &
+           sun_dir_lon,                                                                            &
+           startwall,                                                                              &
+           startland,                                                                              &
+           skyvf,                                                                                  &
+           skyvft,                                                                                 &
+           time_radiation,                                                                         &
+           unscheduled_radiation_calls
 …
 !------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> This subroutine controls the calls of the radiation schemes
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_control
+       IMPLICIT NONE
+       IF ( debug_output_timestep )  CALL debug_message( 'radiation_control', 'start' )
+       SELECT CASE ( TRIM( radiation_scheme ) )
+          CASE ( 'constant' )
+             CALL radiation_constant
+          CASE ( 'clear-sky' )
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE radiation_control
+    IMPLICIT NONE
+    IF ( debug_output_timestep )  CALL debug_message( 'radiation_control', 'start' )
+    SELECT CASE ( TRIM( radiation_scheme ) )
+       CASE ( 'constant' )
+          CALL radiation_constant
+       CASE ( 'clear-sky' )
+          CALL radiation_clearsky
+       CASE ( 'rrtmg' )
+          CALL radiation_rrtmg
+       CASE ( 'external' )
+!
+!--       During spinup apply clear-sky model
+          IF ( time_since_reference_point < 0.0_wp )  THEN
              CALL radiation_clearsky
+          CASE ( 'rrtmg' )
+             CALL radiation_rrtmg
+          CASE ( 'external' )
+!
+!--          During spinup apply clear-sky model
+             IF ( time_since_reference_point < 0.0_wp )  THEN
+                CALL radiation_clearsky
+             ELSE
+                CALL radiation_external
+             ENDIF
+          CASE DEFAULT
+       END SELECT
+       IF ( debug_output_timestep )  CALL debug_message( 'radiation_control', 'end' )
+    END SUBROUTINE radiation_control
+!------------------------------------------------------------------------------!
+          ELSE
+             CALL radiation_external
+          ENDIF
+       CASE DEFAULT
+    END SELECT
+    IF ( debug_output_timestep )  CALL debug_message( 'radiation_control', 'end' )
+ END SUBROUTINE radiation_control
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Check data output for radiation model
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_check_data_output( variable, unit, i, ilen, k )
+       USE control_parameters,                                                 &
+           ONLY: data_output, message_string
+       IMPLICIT NONE
+       LOGICAL                      ::  directional
+       CHARACTER(LEN=*)             ::  unit        !<
+       CHARACTER(LEN=*)             ::  variable    !<
+       CHARACTER(LEN=varnamelength) ::  var         !< TRIM(variable)
+       INTEGER(iwp)                 ::  i, k
+       INTEGER(iwp)                 ::  ilast_word
+       INTEGER(iwp)                 ::  ilen
+       var = TRIM(variable)
+!
+!--    Identify directional variables
+       ilast_word = SCAN(var, '_', back=.TRUE.)
+       IF ( ilast_word > 0 )  THEN
+          SELECT CASE ( var(ilast_word:) )
+             CASE ( '_roof', '_south', '_north', '_west', '_east' )
+                directional = .TRUE.
+                write(9, *) 'vardir', var
+                flush(9)
+                var = var(1:ilast_word-1)
+             CASE DEFAULT
+                directional = .FALSE.
+                write(9, *) 'varnd', var
+                flush(9)
+          END SELECT
+       ELSE
+          directional = .FALSE.
+       END IF
+       IF ( directional )  THEN
+          SELECT CASE ( var )
+             CASE ( 'rtm_rad_net', 'rtm_rad_insw', 'rtm_rad_inlw',             &
+                    'rtm_rad_inswdir', 'rtm_rad_inswdif', 'rtm_rad_inswref',   &
+                    'rtm_rad_inlwdif', 'rtm_rad_inlwref', 'rtm_rad_outsw',     &
+                    'rtm_rad_outlw', 'rtm_rad_ressw', 'rtm_rad_reslw' )
+                IF ( .NOT.  radiation ) THEN
+                   message_string = 'output of "' // var // '" require'        &
+                                    // 's radiation = .TRUE.'
+                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'W/m2'
+             CASE ( 'rtm_svf', 'rtm_dif', 'rtm_skyvf', 'rtm_skyvft',           &
+                    'rtm_surfalb', 'rtm_surfemis' )
+                IF ( .NOT.  radiation ) THEN
+                   message_string = 'output of "' // var // '" require'&
+                                    // 's radiation = .TRUE.'
+                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = '1'
+             CASE DEFAULT
+                unit = 'illegal'
+          END SELECT
+       ELSE
+          SELECT CASE ( var )
+             CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out',     &
+                    'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out'  )
+                IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
+                   message_string = '"output of "' // var // '" requi' //       &
+                                    'res radiation = .TRUE. and ' //            &
+                                    'radiation_scheme = "rrtmg"'
+                   CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'K/h'
+             CASE ( 'rad_net*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',      &
+                    'rad_sw_out*' )
+                IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
+                   message_string = 'illegal value for data_output: "' //       &
+                                    var // '" & only 2d-horizontal ' //         &
+                                    'cross sections are allowed for this value'
+                   CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'W/m2'
+             CASE ( 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*' )
+                IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
+                   message_string = 'illegal value for data_output: "' //       &
+                                    var // '" & only 2d-horizontal ' //         &
+                                    'cross sections are allowed for this value'
+                   CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
+                ENDIF
+                IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
+                   message_string = 'output of "' // var // '" require'         &
+                                    // 's radiation = .TRUE. and radiation_sch' &
+                                    // 'eme = "rrtmg"'
+                   CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = ''
+             CASE ( 'rtm_rad_pc_inlw', 'rtm_rad_pc_insw', 'rtm_rad_pc_inswdir', &
+                    'rtm_rad_pc_inswdif', 'rtm_rad_pc_inswref')
+                IF ( .NOT.  radiation ) THEN
+                   message_string = 'output of "' // var // '" require'         &
+                                    // 's radiation = .TRUE.'
+                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'W'
+             CASE ( 'rtm_mrt', 'rtm_mrt_sw', 'rtm_mrt_lw'  )
+                IF ( .NOT.  radiation ) THEN
+                   message_string = 'output of "' // var // '" require'         &
+                                    // 's radiation = .TRUE.'
+                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+                ENDIF
+                IF ( mrt_nlevels == 0 ) THEN
+                   message_string = 'output of "' // var // '" require'         &
+                                    // 's mrt_nlevels > 0'
+                   CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
+                ENDIF
+                IF ( var == 'rtm_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
+                   message_string = 'output of "' // var // '" require'         &
+                                    // 's rtm_mrt_sw = .TRUE.'
+                   CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
+                ENDIF
+                IF ( var == 'rtm_mrt' ) THEN
+                   unit = 'K'
+                ELSE
+                   unit = 'W m-2'
+                ENDIF
+             CASE DEFAULT
+                unit = 'illegal'
+          END SELECT
+       END IF
+    END SUBROUTINE radiation_check_data_output
+!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE radiation_check_data_output( variable, unit, i, ilen, k )
+    USE control_parameters,                                                                        &
+        ONLY: data_output,                                                                         &
+              message_string
+    IMPLICIT NONE
+    CHARACTER(LEN=*)             ::  unit         !<
+    CHARACTER(LEN=*)             ::  variable     !<
+    CHARACTER(LEN=varnamelength) ::  var          !< TRIM(variable)
+    INTEGER(iwp)                 ::  i, k         !<
+    INTEGER(iwp)                 ::  ilast_word   !<
+    INTEGER(iwp)                 ::  ilen         !<
+    LOGICAL                      ::  directional  !<
+    var = TRIM( variable )
+!
+!-- Identify directional variables
+    ilast_word = SCAN( var, '_', back = .TRUE. )
+    IF ( ilast_word > 0 )  THEN
+       SELECT CASE ( var(ilast_word:) )
+          CASE ( '_roof', '_south', '_north', '_west', '_east' )
+             directional = .TRUE.
+             WRITE( 9, * ) 'vardir', var
+             FLUSH( 9 )
+             var = var(1:ilast_word-1)
+          CASE DEFAULT
+             directional = .FALSE.
+             WRITE( 9, * ) 'varnd', var
+             FLUSH( 9 )
+       END SELECT
+    ELSE
+       directional = .FALSE.
+    END IF
+    IF ( directional )  THEN
+       SELECT CASE ( var )
+          CASE ( 'rtm_rad_net', 'rtm_rad_insw', 'rtm_rad_inlw', 'rtm_rad_inswdir',                 &
+                 'rtm_rad_inswdif', 'rtm_rad_inswref', 'rtm_rad_inlwdif', 'rtm_rad_inlwref',       &
+                 'rtm_rad_outsw', 'rtm_rad_outlw', 'rtm_rad_ressw', 'rtm_rad_reslw' )
+             IF ( .NOT.  radiation )  THEN
+                message_string = 'output of "' // var // '" requires radiation = .TRUE.'
+                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = 'W/m2'
+          CASE ( 'rtm_svf', 'rtm_dif', 'rtm_skyvf', 'rtm_skyvft', 'rtm_surfalb', 'rtm_surfemis' )
+             IF ( .NOT.  radiation )  THEN
+                message_string = 'output of "' // var // '" requires radiation = .TRUE.'
+                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = '1'
+          CASE DEFAULT
+             unit = 'illegal'
+       END SELECT
+    ELSE
+       SELECT CASE ( var )
+          CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out', 'rad_sw_cs_hr',           &
+                 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out'  )
+             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
+                message_string = '"output of "' // var // '" requires radiation = .TRUE. and ' //  &
+                                 'radiation_scheme = "rrtmg"'
+                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = 'K/h'
+          CASE ( 'rad_net*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*', 'rad_sw_out*' )
+             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
+                message_string = 'illegal value for data_output: "' // var //                      &
+                                 '" & only 2d-horizontal cross sections are allowed for this value'
+                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = 'W/m2'
+          CASE ( 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*' )
+             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
+                message_string = 'illegal value for data_output: "' // var //                      &
+                                 '" & only 2d-horizontal cross sections are allowed for this value'
+                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
+                message_string = 'output of "' // var // '" requires radiation = .TRUE. ' //       &
+                                 'and radiation_scheme = "rrtmg"'
+                CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = ''
+          CASE ( 'rtm_rad_pc_inlw', 'rtm_rad_pc_insw', 'rtm_rad_pc_inswdir', 'rtm_rad_pc_inswdif', &
+                 'rtm_rad_pc_inswref')
+             IF ( .NOT.  radiation )  THEN
+                message_string = 'output of "' // var // '" requires radiation = .TRUE.'
+                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = 'W'
+          CASE ( 'rtm_mrt', 'rtm_mrt_sw', 'rtm_mrt_lw'  )
+             IF ( .NOT.  radiation )  THEN
+                message_string = 'output of "' // var // '" requires radiation = .TRUE.'
+                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( mrt_nlevels == 0 )  THEN
+                message_string = 'output of "' // var // '" requires mrt_nlevels > 0'
+                CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( var == 'rtm_mrt_sw'  .AND.  .NOT. mrt_include_sw )  THEN
+                message_string = 'output of "' // var // '" requires rtm_mrt_sw = .TRUE.'
+                CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( var == 'rtm_mrt' )  THEN
+                unit = 'K'
+             ELSE
+                unit = 'W m-2'
+             ENDIF
+          CASE DEFAULT
+             unit = 'illegal'
+       END SELECT
+    END IF
+ END SUBROUTINE radiation_check_data_output
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Set module-specific timeseries units and labels
 !------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE radiation_check_data_output_ts( dots_max, dots_num )
     INTEGER(iwp),      INTENT(IN)     ::  dots_max
     INTEGER(iwp),      INTENT(INOUT)  ::  dots_num
+    INTEGER(iwp), INTENT(IN)    ::  dots_max  !<
+    INTEGER(iwp), INTENT(INOUT) ::  dots_num  !<
+!
 …
+!
+!-- Temporary solution to add LSM and radiation time series to the default
+!-- output
+!-- Temporary solution to add LSM and radiation time series to the default output
     IF ( land_surface  .OR.  radiation )  THEN
        IF ( TRIM( radiation_scheme ) == 'rrtmg' )  THEN
 …
  END SUBROUTINE radiation_check_data_output_ts
 !------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Check data output of profiles for radiation model
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit,      &
+               dopr_unit )
+       USE arrays_3d,                                                          &
+           ONLY: zu
+       USE control_parameters,                                                 &
+           ONLY: data_output_pr, message_string
+       USE indices
+       USE profil_parameter
+       USE statistics
+       IMPLICIT NONE
+       CHARACTER (LEN=*) ::  unit      !<
+       CHARACTER (LEN=*) ::  variable  !<
+       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
+       INTEGER(iwp) ::  var_count     !<
+       SELECT CASE ( TRIM( variable ) )
+         CASE ( 'rad_net' )
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 99
+                dopr_unit  = 'W/m2'
+                hom(:,2,99,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+          CASE ( 'rad_lw_in' )
+             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 100
+                dopr_unit  = 'W/m2'
+                hom(:,2,100,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+          CASE ( 'rad_lw_out' )
+             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 101
+                dopr_unit  = 'W/m2'
+                hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+          CASE ( 'rad_sw_in' )
+             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 102
+                dopr_unit  = 'W/m2'
+                hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+          CASE ( 'rad_sw_out')
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 103
+                dopr_unit  = 'W/m2'
+                hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+          CASE ( 'rad_lw_cs_hr' )
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 104
+                dopr_unit  = 'K/h'
+                hom(:,2,104,:)  = SPREAD( zu, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+          CASE ( 'rad_lw_hr' )
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 105
+                dopr_unit  = 'K/h'
+                hom(:,2,105,:)  = SPREAD( zu, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+          CASE ( 'rad_sw_cs_hr' )
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 106
+                dopr_unit  = 'K/h'
+                hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+          CASE ( 'rad_sw_hr' )
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 107
+                dopr_unit  = 'K/h'
+                hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+          CASE DEFAULT
+             unit = 'illegal'
+       END SELECT
+    END SUBROUTINE radiation_check_data_output_pr
+!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, dopr_unit )
+    USE arrays_3d,                                                                                 &
+        ONLY: zu
+    USE control_parameters,                                                                        &
+        ONLY: data_output_pr,                                                                      &
+              message_string
+    USE indices
+    USE profil_parameter
+    USE statistics
+    IMPLICIT NONE
+    CHARACTER(LEN=*) ::  unit       !<
+    CHARACTER(LEN=*) ::  variable   !<
+    CHARACTER(LEN=*) ::  dopr_unit  !< local value of dopr_unit
+    INTEGER(iwp) ::  var_count  !<
+    SELECT CASE ( TRIM( variable ) )
+      CASE ( 'rad_net' )
+          IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )  THEN
+             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
+                              'not available for radiation = .FALSE. or ' //                       &
+                              'radiation_scheme = "constant"'
+             CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+          ELSE
+             dopr_index(var_count) = 99
+             dopr_unit  = 'W/m2'
+             hom(:,2,99,:)  = SPREAD( zw, 2, statistic_regions + 1 )
+             unit = dopr_unit
+          ENDIF
+       CASE ( 'rad_lw_in' )
+          IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' )  THEN
+             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
+                              'not available for radiation = .FALSE. or ' //                       &
+                              'radiation_scheme = "constant"'
+             CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+          ELSE
+             dopr_index(var_count) = 100
+             dopr_unit  = 'W/m2'
+             hom(:,2,100,:)  = SPREAD( zw, 2, statistic_regions + 1 )
+             unit = dopr_unit
+          ENDIF
+       CASE ( 'rad_lw_out' )
+          IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' )  THEN
+             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
+                              'not available for radiation = .FALSE. or ' //                       &
+                              'radiation_scheme = "constant"'
+             CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+          ELSE
+             dopr_index(var_count) = 101
+             dopr_unit  = 'W/m2'
+             hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions + 1 )
+             unit = dopr_unit
+          ENDIF
+       CASE ( 'rad_sw_in' )
+          IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' )  THEN
+             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
+                              'not available for radiation = .FALSE. or ' //                       &
+                              'radiation_scheme = "constant"'
+             CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+          ELSE
+             dopr_index(var_count) = 102
+             dopr_unit  = 'W/m2'
+             hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions + 1 )
+             unit = dopr_unit
+          ENDIF
+       CASE ( 'rad_sw_out')
+          IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )  THEN
+             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
+                              'not available for radiation = .FALSE. or ' //                       &
+                              'radiation_scheme = "constant"'
+             CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+          ELSE
+             dopr_index(var_count) = 103
+             dopr_unit  = 'W/m2'
+             hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions + 1 )
+             unit = dopr_unit
+          ENDIF
+       CASE ( 'rad_lw_cs_hr' )
+          IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )  THEN
+             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
+                              'not available for radiation = .FALSE. or ' //                       &
+                              'radiation_scheme /= "rrtmg"'
+             CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+          ELSE
+             dopr_index(var_count) = 104
+             dopr_unit  = 'K/h'
+             hom(:,2,104,:)  = SPREAD( zu, 2, statistic_regions + 1 )
+             unit = dopr_unit
+          ENDIF
+       CASE ( 'rad_lw_hr' )
+          IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )  THEN
+             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
+                              'not available for radiation = .FALSE. or ' //                       &
+                              'radiation_scheme /= "rrtmg"'
+             CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+          ELSE
+             dopr_index(var_count) = 105
+             dopr_unit  = 'K/h'
+             hom(:,2,105,:)  = SPREAD( zu, 2, statistic_regions + 1 )
+             unit = dopr_unit
+          ENDIF
+       CASE ( 'rad_sw_cs_hr' )
+          IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )  THEN
+             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
+                              'not available for radiation = .FALSE. or ' //                       &
+                              'radiation_scheme /= "rrtmg"'
+             CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+          ELSE
+             dopr_index(var_count) = 106
+             dopr_unit  = 'K/h'
+             hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions + 1 )
+             unit = dopr_unit
+          ENDIF
+       CASE ( 'rad_sw_hr' )
+          IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )  THEN
+             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
+                              'not available for radiation = .FALSE. or ' //                       &
+                              'radiation_scheme /= "rrtmg"'
+             CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+          ELSE
+             dopr_index(var_count) = 107
+             dopr_unit  = 'K/h'
+             hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions + 1 )
+             unit = dopr_unit
+          ENDIF
+       CASE DEFAULT
+          unit = 'illegal'
+    END SELECT
+ END SUBROUTINE radiation_check_data_output_pr
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Check parameters routine for radiation model
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_check_parameters
+       USE control_parameters,                                                 &
+           ONLY: land_surface, message_string, urban_surface
+       USE netcdf_data_input_mod,                                              &
+           ONLY:  input_pids_static
+       IMPLICIT NONE
+!
+!--    In case no urban-surface or land-surface model is applied, usage of
+!--    a radiation model make no sense.
+       IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
+          message_string = 'Usage of radiation module is only allowed if ' //  &
+                           'land-surface and/or urban-surface model is applied.'
+          CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
+       ENDIF
+       IF ( radiation_scheme /= 'constant'   .AND.                             &
+            radiation_scheme /= 'clear-sky'  .AND.                             &
+            radiation_scheme /= 'rrtmg'      .AND.                             &
+            radiation_scheme /= 'external' )  THEN
+          message_string = 'unknown radiation_scheme = '//                     &
+                           TRIM( radiation_scheme )
+          CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
+       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE radiation_check_parameters
+    USE control_parameters,                                                                        &
+        ONLY: land_surface,                                                                        &
+              message_string,                                                                      &
+              urban_surface
+    USE netcdf_data_input_mod,                                                                     &
+        ONLY:  input_pids_static
+    IMPLICIT NONE
+!
+!-- In case no urban-surface or land-surface model is applied, usage of a radiation model makes no
+!-- sense.
+    IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
+       message_string = 'Usage of radiation module is only allowed if ' //                         &
+                        'land-surface and/or urban-surface model is applied.'
+       CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
+    ENDIF
+    IF ( radiation_scheme /= 'constant'  .AND.  radiation_scheme /= 'clear-sky'  .AND.             &
+         radiation_scheme /= 'rrtmg'  .AND.  radiation_scheme /= 'external' )  THEN
+       message_string = 'unknown radiation_scheme = '// TRIM( radiation_scheme )
+       CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
+    ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
 #if ! defined ( __rrtmg )
+          message_string = 'radiation_scheme = "rrtmg" requires ' //           &
+                           'compilation of PALM with pre-processor ' //        &
+                           'directive -D__rrtmg'
+          CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
+       message_string = 'radiation_scheme = "rrtmg" requires compilation of PALM with ' //         &
+                        'pre-processor directive -D__rrtmg'
+       CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
 #endif
 #if defined ( __rrtmg ) && ! defined( __netcdf )
+          message_string = 'radiation_scheme = "rrtmg" requires ' //           &
+                           'the use of NetCDF (preprocessor directive ' //     &
+                           '-D__netcdf'
+          CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
+       message_string = 'radiation_scheme = "rrtmg" requires the use of NetCDF ' //                &
+                        '(preprocessor directive -D__netcdf'
+       CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
 #endif
+       ENDIF
+!
+!--    Checks performed only if data is given via namelist only.
+       IF ( .NOT. input_pids_static )  THEN
+          IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.          &
+               radiation_scheme == 'clear-sky')  THEN
+             message_string = 'radiation_scheme = "clear-sky" in combination'//&
+                              'with albedo_type = 0 requires setting of'//     &
+                              'albedo /= 9999999.9'
+             CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.     &
+             ( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
+          .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp&
+             ) ) THEN
+             message_string = 'radiation_scheme = "rrtmg" in combination' //   &
+                              'with albedo_type = 0 requires setting of ' //   &
+                              'albedo_lw_dif /= 9999999.9' //                  &
+                              'albedo_lw_dir /= 9999999.9' //                  &
+                              'albedo_sw_dif /= 9999999.9 and' //              &
+                              'albedo_sw_dir /= 9999999.9'
+             CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
+          ENDIF
+       ENDIF
+!
+!--    Parallel rad_angular_discretization without raytrace_mpi_rma is not implemented
+!--    Serial mode does not allow mpi_rma
+    ENDIF
+!
+!-- Checks performed if data is given via namelist only.
+    IF ( .NOT. input_pids_static )  THEN
+       IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.                                 &
+            radiation_scheme == 'clear-sky')  THEN
+          message_string = 'radiation_scheme = "clear-sky" in combination with albedo_type = 0 ' //&
+                           'requires setting of albedo /= 9999999.9'
+          CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
+       ENDIF
+       IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.                            &
+          ( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp                       &
+       .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp ) )  THEN
+          message_string = 'radiation_scheme = "rrtmg" in combination with albedo_type = 0 ' //    &
+                           'requires setting of albedo_lw_dif /= 9999999.9' //                     &
+                           'albedo_lw_dir /= 9999999.9' //                                         &
+                           'albedo_sw_dif /= 9999999.9 and albedo_sw_dir /= 9999999.9'
+          CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
+       ENDIF
+    ENDIF
+!
+!-- Parallel rad_angular_discretization without raytrace_mpi_rma is not implemented. Serial mode does
+!-- not allow mpi_rma
 #if defined( __parallel )
+       IF ( rad_angular_discretization  .AND.  .NOT. raytrace_mpi_rma )  THEN
+          message_string = 'rad_angular_discretization can only be used ' //  &
+                           'together with raytrace_mpi_rma or when ' //  &
+                           'no parallelization is applied.'
+          CALL message( 'readiation_check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
+       ENDIF
+    IF ( rad_angular_discretization  .AND.  .NOT. raytrace_mpi_rma )  THEN
+       message_string = 'rad_angular_discretization can only be used together with ' //            &
+                        'raytrace_mpi_rma or when no parallelization is applied.'
+       CALL message( 'readiation_check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
+    ENDIF
 #else
        IF ( raytrace_mpi_rma )  THEN
           message_string = 'raytrace_mpi_rma = .T. not allowed in serial mode'
           CALL message( 'readiation_check_parameters', 'PA0710', 1, 2, 0, 6, 0 )
        ENDIF
+    IF ( raytrace_mpi_rma )  THEN
+       message_string = 'raytrace_mpi_rma = .T. not allowed in serial mode'
+       CALL message( 'readiation_check_parameters', 'PA0710', 1, 2, 0, 6, 0 )
+    ENDIF
 #endif
+       IF ( cloud_droplets  .AND.   radiation_scheme == 'rrtmg'  .AND.         &
+            average_radiation ) THEN
+          message_string = 'average_radiation = .T. with radiation_scheme'//   &
+                           '= "rrtmg" in combination cloud_droplets = .T.'//   &
+                           'is not implementd'
+          CALL message( 'check_parameters', 'PA0560', 1, 2, 0, 6, 0 )
+       ENDIF
+!
+!--    Incialize svf normalization reporting histogram
+       svfnorm_report_num = 1
+       DO WHILE ( svfnorm_report_thresh(svfnorm_report_num) < 1e20_wp          &
+                   .AND. svfnorm_report_num <= 30 )
+          svfnorm_report_num = svfnorm_report_num + 1
+       ENDDO
+       svfnorm_report_num = svfnorm_report_num - 1
+!
+!--    Check for dt_radiation
+       IF ( dt_radiation <= 0.0 )  THEN
+          message_string = 'dt_radiation must be > 0.0'
+          CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 )
+       ENDIF
+    END SUBROUTINE radiation_check_parameters
+!------------------------------------------------------------------------------!
+    IF ( cloud_droplets  .AND.   radiation_scheme == 'rrtmg'  .AND.  average_radiation )  THEN
+       message_string = 'average_radiation = .T. with radiation_scheme = "rrtmg" ' //              &
+                        'in combination cloud_droplets = .T. is not implementd'
+       CALL message( 'check_parameters', 'PA0560', 1, 2, 0, 6, 0 )
+    ENDIF
+!
+!-- Initialize svf normalization reporting histogram
+    svfnorm_report_num = 1
+    DO WHILE ( svfnorm_report_thresh(svfnorm_report_num) < 1E20_wp .AND. svfnorm_report_num <= 30 )
+       svfnorm_report_num = svfnorm_report_num + 1
+    ENDDO
+    svfnorm_report_num = svfnorm_report_num - 1
+!
+!-- Check for dt_radiation
+    IF ( dt_radiation <= 0.0 )  THEN
+       message_string = 'dt_radiation must be > 0.0'
+       CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 )
+    ENDIF
+ END SUBROUTINE radiation_check_parameters
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Initialization of the radiation model and Radiative Transfer Model
 !------------------------------------------------------------------------------!
     SUBROUTINE radiation_init
        IMPLICIT NONE
        INTEGER(iwp) ::  i         !< running index x-direction
        INTEGER(iwp) ::  is        !< running index for input surface elements
        INTEGER(iwp) ::  ioff      !< offset in x between surface element reference grid point in atmosphere and actual surface
        INTEGER(iwp) ::  j         !< running index y-direction
        INTEGER(iwp) ::  joff      !< offset in y between surface element reference grid point in atmosphere and actual surface
        INTEGER(iwp) ::  k         !< running index z-direction
        INTEGER(iwp) ::  l         !< running index for orientation of vertical surfaces
        INTEGER(iwp) ::  m         !< running index for surface elements
        INTEGER(iwp) ::  ntime = 0 !< number of available external radiation timesteps
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE radiation_init
+    IMPLICIT NONE
+    INTEGER(iwp) ::  i          !< running index x-direction
+    INTEGER(iwp) ::  is         !< running index for input surface elements
+    INTEGER(iwp) ::  ioff       !< offset in x between surface element reference grid point in atmosphere and actual surface
+    INTEGER(iwp) ::  j          !< running index y-direction
+    INTEGER(iwp) ::  joff       !< offset in y between surface element reference grid point in atmosphere and actual surface
+    INTEGER(iwp) ::  k          !< running index z-direction
+    INTEGER(iwp) ::  l          !< running index for orientation of vertical surfaces
+    INTEGER(iwp) ::  m          !< running index for surface elements
+    INTEGER(iwp) ::  ntime = 0  !< number of available external radiation timesteps
 #if defined( __rrtmg )
        INTEGER(iwp) ::  ind_type  !< running index for subgrid-surface tiles
+    INTEGER(iwp) :: ind_type   !< running index for subgrid-surface tiles
 #endif
+       LOGICAL      ::  radiation_input_root_domain !< flag indicating the existence of a dynamic input file for the root domain
+       IF ( debug_output )  CALL debug_message( 'radiation_init', 'start' )
+!
+!--    Activate radiation_interactions according to the existence of vertical surfaces and/or trees
+!      or if biometeorology output is required for flat surfaces.
+!--    The namelist parameter radiation_interactions_on can override this behavior.
+!--    (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in
+!--    init_surface_arrays.)
+       IF ( radiation_interactions_on )  THEN
+          IF ( vertical_surfaces_exist  .OR.  plant_canopy  .OR.  biometeorology )  THEN
+             radiation_interactions    = .TRUE.
+             average_radiation         = .TRUE.
+          ELSE
+             radiation_interactions_on = .FALSE.   !< reset namelist parameter: no interactions
+                                                   !< calculations necessary in case of flat surface
+    LOGICAL ::  radiation_input_root_domain  !< flag indicating the existence of a dynamic input file for the root domain
+    IF ( debug_output )  CALL debug_message( 'radiation_init', 'start' )
+!
+!-- Activate radiation_interactions according to the existence of vertical surfaces and/or trees
+!   or if biometeorology output is required for flat surfaces.
+!-- The namelist parameter radiation_interactions_on can override this behavior. (This check cannot
+!-- be performed in check_parameters, because vertical_surfaces_exist is first set in
+!-- init_surface_arrays.)
+    IF ( radiation_interactions_on )  THEN
+       IF ( vertical_surfaces_exist  .OR.  plant_canopy  .OR.  biometeorology )  THEN
+          radiation_interactions    = .TRUE.
+          average_radiation         = .TRUE.
+       ELSE
+          radiation_interactions_on = .FALSE.  !< reset namelist parameter: no interactions
+                                               !< calculations necessary in case of flat surface
+       ENDIF
+    ELSEIF ( vertical_surfaces_exist  .OR.  plant_canopy  .OR.  biometeorology )  THEN
+       message_string = 'radiation_interactions_on is set to .FALSE. although vertical ' //        &
+                        'surfaces and/or trees or biometeorology exist is ON. The model will ' //  &
+                        'run without RTM (no shadows, no radiation reflections)'
+       CALL message( 'init_3d_model', 'PA0348', 0, 1, 0, 6, 0 )
+    ENDIF
+!
+!-- Precalculate some time constants
+    d_hours_day    = 1.0_wp / REAL( hours_per_day, KIND = wp )
+    d_seconds_hour = 1.0_wp / seconds_per_hour
+!
+!-- If required, initialize radiation interactions between surfaces via sky-view factors. This must
+!-- be done before radiation is initialized.
+    IF ( radiation_interactions )  CALL radiation_interaction_init
+!
+!-- Allocate array for storing the surface net radiation
+    IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_net )  .AND.  surf_lsm_h%ns > 0 )  THEN
+       ALLOCATE( surf_lsm_h%rad_net(1:surf_lsm_h%ns) )
+       surf_lsm_h%rad_net = 0.0_wp
+    ENDIF
+    IF ( .NOT. ALLOCATED ( surf_usm_h%rad_net )  .AND.  surf_usm_h%ns > 0 )  THEN
+       ALLOCATE( surf_usm_h%rad_net(1:surf_usm_h%ns) )
+       surf_usm_h%rad_net = 0.0_wp
+    ENDIF
+    DO  l = 0, 3
+       IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_net )  .AND.  surf_lsm_v(l)%ns > 0 )  THEN
+          ALLOCATE( surf_lsm_v(l)%rad_net(1:surf_lsm_v(l)%ns) )
+          surf_lsm_v(l)%rad_net = 0.0_wp
+       ENDIF
+       IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_net )  .AND.  surf_usm_v(l)%ns > 0 )  THEN
+          ALLOCATE( surf_usm_v(l)%rad_net(1:surf_usm_v(l)%ns) )
+          surf_usm_v(l)%rad_net = 0.0_wp
+       ENDIF
+    ENDDO
+!
+!-- Allocate array for storing the surface longwave (out) radiation change
+    IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_lw_out_change_0 )  .AND.  surf_lsm_h%ns > 0 )  THEN
+       ALLOCATE( surf_lsm_h%rad_lw_out_change_0(1:surf_lsm_h%ns) )
+       surf_lsm_h%rad_lw_out_change_0 = 0.0_wp
+    ENDIF
+    IF ( .NOT. ALLOCATED ( surf_usm_h%rad_lw_out_change_0 )  .AND.  surf_usm_h%ns > 0 ) THEN
+       ALLOCATE( surf_usm_h%rad_lw_out_change_0(1:surf_usm_h%ns) )
+       surf_usm_h%rad_lw_out_change_0 = 0.0_wp
+    ENDIF
+    DO  l = 0, 3
+       IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_lw_out_change_0 )  .AND.  surf_lsm_v(l)%ns > 0 )   &
+       THEN
+          ALLOCATE( surf_lsm_v(l)%rad_lw_out_change_0(1:surf_lsm_v(l)%ns) )
+          surf_lsm_v(l)%rad_lw_out_change_0 = 0.0_wp
+       ENDIF
+       IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_lw_out_change_0 )  .AND.  surf_usm_v(l)%ns > 0 )   &
+       THEN
+          ALLOCATE( surf_usm_v(l)%rad_lw_out_change_0(1:surf_usm_v(l)%ns) )
+          surf_usm_v(l)%rad_lw_out_change_0 = 0.0_wp
+       ENDIF
+    ENDDO
+!
+!-- Allocate surface arrays for incoming/outgoing short/longwave radiation
+    IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_sw_in )  .AND.  surf_lsm_h%ns > 0 )  THEN
+       ALLOCATE( surf_lsm_h%rad_sw_in(1:surf_lsm_h%ns)  )
+       ALLOCATE( surf_lsm_h%rad_sw_out(1:surf_lsm_h%ns) )
+       ALLOCATE( surf_lsm_h%rad_sw_dir(1:surf_lsm_h%ns) )
+       ALLOCATE( surf_lsm_h%rad_sw_dif(1:surf_lsm_h%ns) )
+       ALLOCATE( surf_lsm_h%rad_sw_ref(1:surf_lsm_h%ns) )
+       ALLOCATE( surf_lsm_h%rad_sw_res(1:surf_lsm_h%ns) )
+       ALLOCATE( surf_lsm_h%rad_lw_in(1:surf_lsm_h%ns)  )
+       ALLOCATE( surf_lsm_h%rad_lw_out(1:surf_lsm_h%ns) )
+       ALLOCATE( surf_lsm_h%rad_lw_dif(1:surf_lsm_h%ns) )
+       ALLOCATE( surf_lsm_h%rad_lw_ref(1:surf_lsm_h%ns) )
+       ALLOCATE( surf_lsm_h%rad_lw_res(1:surf_lsm_h%ns) )
+       surf_lsm_h%rad_sw_in  = 0.0_wp
+       surf_lsm_h%rad_sw_out = 0.0_wp
+       surf_lsm_h%rad_sw_dir = 0.0_wp
+       surf_lsm_h%rad_sw_dif = 0.0_wp
+       surf_lsm_h%rad_sw_ref = 0.0_wp
+       surf_lsm_h%rad_sw_res = 0.0_wp
+       surf_lsm_h%rad_lw_in  = 0.0_wp
+       surf_lsm_h%rad_lw_out = 0.0_wp
+       surf_lsm_h%rad_lw_dif = 0.0_wp
+       surf_lsm_h%rad_lw_ref = 0.0_wp
+       surf_lsm_h%rad_lw_res = 0.0_wp
+    ENDIF
+    IF ( .NOT. ALLOCATED ( surf_usm_h%rad_sw_in )  .AND.  surf_usm_h%ns > 0 ) THEN
+       ALLOCATE( surf_usm_h%rad_sw_in(1:surf_usm_h%ns)  )
+       ALLOCATE( surf_usm_h%rad_sw_out(1:surf_usm_h%ns) )
+       ALLOCATE( surf_usm_h%rad_sw_dir(1:surf_usm_h%ns) )
+       ALLOCATE( surf_usm_h%rad_sw_dif(1:surf_usm_h%ns) )
+       ALLOCATE( surf_usm_h%rad_sw_ref(1:surf_usm_h%ns) )
+       ALLOCATE( surf_usm_h%rad_sw_res(1:surf_usm_h%ns) )
+       ALLOCATE( surf_usm_h%rad_lw_in(1:surf_usm_h%ns)  )
+       ALLOCATE( surf_usm_h%rad_lw_out(1:surf_usm_h%ns) )
+       ALLOCATE( surf_usm_h%rad_lw_dif(1:surf_usm_h%ns) )
+       ALLOCATE( surf_usm_h%rad_lw_ref(1:surf_usm_h%ns) )
+       ALLOCATE( surf_usm_h%rad_lw_res(1:surf_usm_h%ns) )
+       surf_usm_h%rad_sw_in  = 0.0_wp
+       surf_usm_h%rad_sw_out = 0.0_wp
+       surf_usm_h%rad_sw_dir = 0.0_wp
+       surf_usm_h%rad_sw_dif = 0.0_wp
+       surf_usm_h%rad_sw_ref = 0.0_wp
+       surf_usm_h%rad_sw_res = 0.0_wp
+       surf_usm_h%rad_lw_in  = 0.0_wp
+       surf_usm_h%rad_lw_out = 0.0_wp
+       surf_usm_h%rad_lw_dif = 0.0_wp
+       surf_usm_h%rad_lw_ref = 0.0_wp
+       surf_usm_h%rad_lw_res = 0.0_wp
+    ENDIF
+    DO  l = 0, 3
+       IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_sw_in )  .AND.  surf_lsm_v(l)%ns > 0 )  THEN
+          ALLOCATE( surf_lsm_v(l)%rad_sw_in(1:surf_lsm_v(l)%ns)  )
+          ALLOCATE( surf_lsm_v(l)%rad_sw_out(1:surf_lsm_v(l)%ns) )
+          ALLOCATE( surf_lsm_v(l)%rad_sw_dir(1:surf_lsm_v(l)%ns) )
+          ALLOCATE( surf_lsm_v(l)%rad_sw_dif(1:surf_lsm_v(l)%ns) )
+          ALLOCATE( surf_lsm_v(l)%rad_sw_ref(1:surf_lsm_v(l)%ns) )
+          ALLOCATE( surf_lsm_v(l)%rad_sw_res(1:surf_lsm_v(l)%ns) )
+          ALLOCATE( surf_lsm_v(l)%rad_lw_in(1:surf_lsm_v(l)%ns)  )
+          ALLOCATE( surf_lsm_v(l)%rad_lw_out(1:surf_lsm_v(l)%ns) )
+          ALLOCATE( surf_lsm_v(l)%rad_lw_dif(1:surf_lsm_v(l)%ns) )
+          ALLOCATE( surf_lsm_v(l)%rad_lw_ref(1:surf_lsm_v(l)%ns) )
+          ALLOCATE( surf_lsm_v(l)%rad_lw_res(1:surf_lsm_v(l)%ns) )
+          surf_lsm_v(l)%rad_sw_in  = 0.0_wp
+          surf_lsm_v(l)%rad_sw_out = 0.0_wp
+          surf_lsm_v(l)%rad_sw_dir = 0.0_wp
+          surf_lsm_v(l)%rad_sw_dif = 0.0_wp
+          surf_lsm_v(l)%rad_sw_ref = 0.0_wp
+          surf_lsm_v(l)%rad_sw_res = 0.0_wp
+          surf_lsm_v(l)%rad_lw_in  = 0.0_wp
+          surf_lsm_v(l)%rad_lw_out = 0.0_wp
+          surf_lsm_v(l)%rad_lw_dif = 0.0_wp
+          surf_lsm_v(l)%rad_lw_ref = 0.0_wp
+          surf_lsm_v(l)%rad_lw_res = 0.0_wp
+       ENDIF
+       IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_sw_in )  .AND.  surf_usm_v(l)%ns > 0 )  THEN
+          ALLOCATE( surf_usm_v(l)%rad_sw_in(1:surf_usm_v(l)%ns)  )
+          ALLOCATE( surf_usm_v(l)%rad_sw_out(1:surf_usm_v(l)%ns) )
+          ALLOCATE( surf_usm_v(l)%rad_sw_dir(1:surf_usm_v(l)%ns) )
+          ALLOCATE( surf_usm_v(l)%rad_sw_dif(1:surf_usm_v(l)%ns) )
+          ALLOCATE( surf_usm_v(l)%rad_sw_ref(1:surf_usm_v(l)%ns) )
+          ALLOCATE( surf_usm_v(l)%rad_sw_res(1:surf_usm_v(l)%ns) )
+          ALLOCATE( surf_usm_v(l)%rad_lw_in(1:surf_usm_v(l)%ns)  )
+          ALLOCATE( surf_usm_v(l)%rad_lw_out(1:surf_usm_v(l)%ns) )
+          ALLOCATE( surf_usm_v(l)%rad_lw_dif(1:surf_usm_v(l)%ns) )
+          ALLOCATE( surf_usm_v(l)%rad_lw_ref(1:surf_usm_v(l)%ns) )
+          ALLOCATE( surf_usm_v(l)%rad_lw_res(1:surf_usm_v(l)%ns) )
+          surf_usm_v(l)%rad_sw_in  = 0.0_wp
+          surf_usm_v(l)%rad_sw_out = 0.0_wp
+          surf_usm_v(l)%rad_sw_dir = 0.0_wp
+          surf_usm_v(l)%rad_sw_dif = 0.0_wp
+          surf_usm_v(l)%rad_sw_ref = 0.0_wp
+          surf_usm_v(l)%rad_sw_res = 0.0_wp
+          surf_usm_v(l)%rad_lw_in  = 0.0_wp
+          surf_usm_v(l)%rad_lw_out = 0.0_wp
+          surf_usm_v(l)%rad_lw_dif = 0.0_wp
+          surf_usm_v(l)%rad_lw_ref = 0.0_wp
+          surf_usm_v(l)%rad_lw_res = 0.0_wp
+       ENDIF
+    ENDDO
+!
+!-- Fix net radiation in case of radiation_scheme = 'constant'
+    IF ( radiation_scheme == 'constant' )  THEN
+       IF ( ALLOCATED( surf_lsm_h%rad_net ) )  surf_lsm_h%rad_net    = net_radiation
+       IF ( ALLOCATED( surf_usm_h%rad_net ) )  surf_usm_h%rad_net    = net_radiation
+!
+!--    Todo: weight with inclination angle
+       DO  l = 0, 3
+          IF ( ALLOCATED( surf_lsm_v(l)%rad_net ) )  surf_lsm_v(l)%rad_net = net_radiation
+          IF ( ALLOCATED( surf_usm_v(l)%rad_net ) )  surf_usm_v(l)%rad_net = net_radiation
+       ENDDO
+!       radiation = .FALSE.
+!
+!-- Calculate orbital constants
+    ELSE
+       decl_1 = SIN( 23.45_wp * pi / 180.0_wp )
+       decl_2 = 2.0_wp * pi / 365.0_wp
+       decl_3 = decl_2 * 81.0_wp
+       lat    = latitude * pi / 180.0_wp
+       lon    = longitude * pi / 180.0_wp
+    ENDIF
+    IF ( radiation_scheme == 'clear-sky'  .OR.                                                     &
+         radiation_scheme == 'constant'   .OR.                                                     &
+         radiation_scheme == 'external' )  THEN
+!
+!--    Allocate arrays for incoming/outgoing short/longwave radiation
+       IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
+          ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
+          rad_sw_in = 0.0_wp
+       ENDIF
+       IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
+          ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
+          rad_sw_out = 0.0_wp
+       ENDIF
+       IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
+          ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
+          rad_lw_in = 0.0_wp
+       ENDIF
+       IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
+          ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
+          rad_lw_out = 0.0_wp
+       ENDIF
+!
+!--    Allocate average arrays for incoming/outgoing short/longwave radiation
+       IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
+          ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
+       ENDIF
+       IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
+          ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
+       ENDIF
+       IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
+          ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
+       ENDIF
+       IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
+          ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
+       ENDIF
+!
+!--    Allocate arrays for broadband albedo, and level 1 initialization via namelist paramter,
+!--    unless not already allocated.
+       IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )  THEN
+          ALLOCATE( surf_lsm_h%albedo(1:surf_lsm_h%ns,0:2) )
+          surf_lsm_h%albedo    = albedo
+       ENDIF
+       IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )  THEN
+          ALLOCATE( surf_usm_h%albedo(1:surf_usm_h%ns,0:2) )
+          surf_usm_h%albedo    = albedo
+       ENDIF
+       DO  l = 0, 3
+          IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )  THEN
+             ALLOCATE( surf_lsm_v(l)%albedo(1:surf_lsm_v(l)%ns,0:2) )
+             surf_lsm_v(l)%albedo = albedo
           ENDIF
+       ELSEIF ( vertical_surfaces_exist  .OR.  plant_canopy  .OR.  biometeorology )  THEN
+          message_string = 'radiation_interactions_on is set to .FALSE. although '     // &
+                           'vertical surfaces and/or trees or biometeorology exist '   // &
+                           'is ON. The model will run without RTM (no shadows, no '    // &
+                           'radiation reflections)'
+          CALL message( 'init_3d_model', 'PA0348', 0, 1, 0, 6, 0 )
+       ENDIF
+!
+!--    Precalculate some time constants
+       d_hours_day    = 1.0_wp / REAL( hours_per_day, KIND = wp )
+       d_seconds_hour = 1.0_wp / seconds_per_hour
+!
+!--    If required, initialize radiation interactions between surfaces
+!--    via sky-view factors. This must be done before radiation is initialized.
+       IF ( radiation_interactions )  CALL radiation_interaction_init
+!
+!--    Allocate array for storing the surface net radiation
+       IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_net )  .AND.                      &
+                  surf_lsm_h%ns > 0  )   THEN
+          ALLOCATE( surf_lsm_h%rad_net(1:surf_lsm_h%ns) )
+          surf_lsm_h%rad_net = 0.0_wp
+       ENDIF
+       IF ( .NOT. ALLOCATED ( surf_usm_h%rad_net )  .AND.                      &
+                  surf_usm_h%ns > 0  )  THEN
+          ALLOCATE( surf_usm_h%rad_net(1:surf_usm_h%ns) )
+          surf_usm_h%rad_net = 0.0_wp
+       ENDIF
+       DO  l = 0, 3
+          IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_net )  .AND.                &
+                     surf_lsm_v(l)%ns > 0  )  THEN
+             ALLOCATE( surf_lsm_v(l)%rad_net(1:surf_lsm_v(l)%ns) )
+             surf_lsm_v(l)%rad_net = 0.0_wp
+          ENDIF
+          IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_net )  .AND.                &
+                     surf_usm_v(l)%ns > 0  )  THEN
+             ALLOCATE( surf_usm_v(l)%rad_net(1:surf_usm_v(l)%ns) )
+             surf_usm_v(l)%rad_net = 0.0_wp
+          IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )  THEN
+             ALLOCATE( surf_usm_v(l)%albedo(1:surf_usm_v(l)%ns,0:2) )
+             surf_usm_v(l)%albedo = albedo
           ENDIF
        ENDDO
+!
+!--    Allocate array for storing the surface longwave (out) radiation change
+       IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_lw_out_change_0 )  .AND.          &
+                  surf_lsm_h%ns > 0  )   THEN
+          ALLOCATE( surf_lsm_h%rad_lw_out_change_0(1:surf_lsm_h%ns) )
+          surf_lsm_h%rad_lw_out_change_0 = 0.0_wp
+       ENDIF
+       IF ( .NOT. ALLOCATED ( surf_usm_h%rad_lw_out_change_0 )  .AND.          &
+                  surf_usm_h%ns > 0  )  THEN
+          ALLOCATE( surf_usm_h%rad_lw_out_change_0(1:surf_usm_h%ns) )
+          surf_usm_h%rad_lw_out_change_0 = 0.0_wp
+       ENDIF
+!
+!--    Level 2 initialization of broadband albedo via given albedo_type.
+!--    Only if albedo_type is non-zero. In case of urban surface and input data is read from ASCII
+!--    file, albedo_type will be zero, so that albedo won't be overwritten.
+       DO  m = 1, surf_lsm_h%ns
+          IF ( surf_lsm_h%albedo_type(m,ind_veg_wall) /= 0 )                                       &
+             surf_lsm_h%albedo(m,ind_veg_wall) =                                                   &
+             albedo_pars(0,surf_lsm_h%albedo_type(m,ind_veg_wall))
+          IF ( surf_lsm_h%albedo_type(m,ind_pav_green) /= 0 )                                      &
+             surf_lsm_h%albedo(m,ind_pav_green) =                                                  &
+                        albedo_pars(0,surf_lsm_h%albedo_type(m,ind_pav_green))
+          IF ( surf_lsm_h%albedo_type(m,ind_wat_win) /= 0 )                                        &
+             surf_lsm_h%albedo(m,ind_wat_win) =                                                    &
+                        albedo_pars(0,surf_lsm_h%albedo_type(m,ind_wat_win))
+       ENDDO
+       DO  m = 1, surf_usm_h%ns
+          IF ( surf_usm_h%albedo_type(m,ind_veg_wall) /= 0 )                                       &
+             surf_usm_h%albedo(m,ind_veg_wall) =                                                   &
+                        albedo_pars(0,surf_usm_h%albedo_type(m,ind_veg_wall))
+          IF ( surf_usm_h%albedo_type(m,ind_pav_green) /= 0 )                                      &
+             surf_usm_h%albedo(m,ind_pav_green) =                                                  &
+                        albedo_pars(0,surf_usm_h%albedo_type(m,ind_pav_green))
+          IF ( surf_usm_h%albedo_type(m,ind_wat_win) /= 0 )                                        &
+             surf_usm_h%albedo(m,ind_wat_win) =                                                    &
+                        albedo_pars(0,surf_usm_h%albedo_type(m,ind_wat_win))
+       ENDDO
        DO  l = 0, 3
+          IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_lw_out_change_0 )  .AND.    &
+                     surf_lsm_v(l)%ns > 0  )  THEN
+             ALLOCATE( surf_lsm_v(l)%rad_lw_out_change_0(1:surf_lsm_v(l)%ns) )
+             surf_lsm_v(l)%rad_lw_out_change_0 = 0.0_wp
+          DO  m = 1, surf_lsm_v(l)%ns
+             IF ( surf_lsm_v(l)%albedo_type(m,ind_veg_wall) /= 0 )                                 &
+                surf_lsm_v(l)%albedo(m,ind_veg_wall) =                                             &
+                     albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_veg_wall))
+             IF ( surf_lsm_v(l)%albedo_type(m,ind_pav_green) /= 0 )                                &
+                surf_lsm_v(l)%albedo(m,ind_pav_green) =                                            &
+                     albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_pav_green))
+             IF ( surf_lsm_v(l)%albedo_type(m,ind_wat_win) /= 0 )                                  &
+                surf_lsm_v(l)%albedo(m,ind_wat_win) =                                              &
+                     albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_wat_win))
+          ENDDO
+          DO  m = 1, surf_usm_v(l)%ns
+             IF ( surf_usm_v(l)%albedo_type(m,ind_veg_wall) /= 0 )                                 &
+                surf_usm_v(l)%albedo(m,ind_veg_wall) =                                             &
+                     albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_veg_wall))
+             IF ( surf_usm_v(l)%albedo_type(m,ind_pav_green) /= 0 )                                &
+                surf_usm_v(l)%albedo(m,ind_pav_green) =                                            &
+                     albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_pav_green))
+             IF ( surf_usm_v(l)%albedo_type(m,ind_wat_win) /= 0 )                                  &
+                surf_usm_v(l)%albedo(m,ind_wat_win) =                                              &
+                     albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_wat_win))
+          ENDDO
+       ENDDO
+!
+!--    Level 3 initialization at grid points where albedo type is zero.
+!--    This case, albedo is taken from file. In case of constant radiation or clear sky, only
+!--    broadband albedo is given.
+       IF ( albedo_pars_f%from_file )  THEN
+!
+!--       Horizontal surfaces
+          DO  m = 1, surf_lsm_h%ns
+             i = surf_lsm_h%i(m)
+             j = surf_lsm_h%j(m)
+             IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
+                surf_lsm_h%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
+                surf_lsm_h%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
+                surf_lsm_h%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
+             ENDIF
+          ENDDO
+          DO  m = 1, surf_usm_h%ns
+             i = surf_usm_h%i(m)
+             j = surf_usm_h%j(m)
+             IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
+                surf_usm_h%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
+                surf_usm_h%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
+                surf_usm_h%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
+             ENDIF
+          ENDDO
+!
+!--       Vertical surfaces
+          DO  l = 0, 3
+             ioff = surf_lsm_v(l)%ioff
+             joff = surf_lsm_v(l)%joff
+             DO  m = 1, surf_lsm_v(l)%ns
+                i = surf_lsm_v(l)%i(m) + ioff
+                j = surf_lsm_v(l)%j(m) + joff
+                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
+                   surf_lsm_v(l)%albedo(m,ind_veg_wall) = albedo_pars_f%pars_xy(0,j,i)
+                   surf_lsm_v(l)%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
+                   surf_lsm_v(l)%albedo(m,ind_wat_win) = albedo_pars_f%pars_xy(0,j,i)
+                ENDIF
+             ENDDO
+             ioff = surf_usm_v(l)%ioff
+             joff = surf_usm_v(l)%joff
+             DO  m = 1, surf_usm_v(l)%ns
+                i = surf_usm_v(l)%i(m) + ioff
+                j = surf_usm_v(l)%j(m) + joff
+                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
+                   surf_usm_v(l)%albedo(m,ind_veg_wall) = albedo_pars_f%pars_xy(0,j,i)
+                   surf_usm_v(l)%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
+                   surf_usm_v(l)%albedo(m,ind_wat_win) = albedo_pars_f%pars_xy(0,j,i)
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDIF
+!
+!--    Read explicit albedo values from building surface pars. If present, they override all less
+!--    specific albedo values and force a albedo_type to zero in order to take effect.
+       IF ( building_surface_pars_f%from_file )  THEN
+          DO  m = 1, surf_usm_h%ns
+             i = surf_usm_h%i(m)
+             j = surf_usm_h%j(m)
+             k = surf_usm_h%k(m)
+!
+!--          Iterate over surfaces in column, check height and orientation
+             DO  is = building_surface_pars_f%index_ji(1,j,i),                                     &
+                      building_surface_pars_f%index_ji(2,j,i)
+                IF ( building_surface_pars_f%coords(4,is) == -surf_usm_h%koff  .AND.               &
+                     building_surface_pars_f%coords(1,is) == k )  THEN
+                   IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=                       &
+                        building_surface_pars_f%fill )  THEN
+                      surf_usm_h%albedo(m,ind_veg_wall) =                                          &
+                               building_surface_pars_f%pars(ind_s_alb_b_wall,is)
+                      surf_usm_h%albedo_type(m,ind_veg_wall) = 0
+                   ENDIF
+                   IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=                        &
+                        building_surface_pars_f%fill )  THEN
+                      surf_usm_h%albedo(m,ind_wat_win) =                                           &
+                               building_surface_pars_f%pars(ind_s_alb_b_win,is)
+                      surf_usm_h%albedo_type(m,ind_wat_win) = 0
+                   ENDIF
+                   IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=                      &
+                        building_surface_pars_f%fill )  THEN
+                      surf_usm_h%albedo(m,ind_pav_green) =                                         &
+                               building_surface_pars_f%pars(ind_s_alb_b_green,is)
+                      surf_usm_h%albedo_type(m,ind_pav_green) = 0
+                   ENDIF
+                   EXIT ! Surface was found and processed
+                ENDIF
+             ENDDO
+          ENDDO
+          DO  l = 0, 3
+             DO  m = 1, surf_usm_v(l)%ns
+                i = surf_usm_v(l)%i(m)
+                j = surf_usm_v(l)%j(m)
+                k = surf_usm_v(l)%k(m)
+!
+!--             Iterate over surfaces in column, check height and orientation
+                DO  is = building_surface_pars_f%index_ji(1,j,i),                                  &
+                         building_surface_pars_f%index_ji(2,j,i)
+                   IF ( building_surface_pars_f%coords(5,is) == -surf_usm_v(l)%joff  .AND.         &
+                        building_surface_pars_f%coords(6,is) == -surf_usm_v(l)%ioff  .AND.         &
+                        building_surface_pars_f%coords(1,is) == k )  THEN
+                      IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=                    &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_v(l)%albedo(m,ind_veg_wall) =                                    &
+                                  building_surface_pars_f%pars(ind_s_alb_b_wall,is)
+                         surf_usm_v(l)%albedo_type(m,ind_veg_wall) = 0
+                      ENDIF
+                      IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=                     &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_v(l)%albedo(m,ind_wat_win) =                                     &
+                                  building_surface_pars_f%pars(ind_s_alb_b_win,is)
+                         surf_usm_v(l)%albedo_type(m,ind_wat_win) = 0
+                      ENDIF
+                      IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=                   &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_v(l)%albedo(m,ind_pav_green) =                                   &
+                                  building_surface_pars_f%pars(ind_s_alb_b_green,is)
+                         surf_usm_v(l)%albedo_type(m,ind_pav_green) = 0
+                      ENDIF
+                      EXIT ! Surface was found and processed
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
+!
+!-- Initialization actions for RRTMG
+    ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
+#if defined ( __rrtmg )
+!
+!--    Allocate albedos for short/longwave radiation, horizontal surfaces for wall/green/window
+!--    (USM) or vegetation/pavement/water surfaces (LSM).
+       ALLOCATE ( surf_lsm_h%aldif(1:surf_lsm_h%ns,0:2)      )
+       ALLOCATE ( surf_lsm_h%aldir(1:surf_lsm_h%ns,0:2)      )
+       ALLOCATE ( surf_lsm_h%asdif(1:surf_lsm_h%ns,0:2)      )
+       ALLOCATE ( surf_lsm_h%asdir(1:surf_lsm_h%ns,0:2)      )
+       ALLOCATE ( surf_lsm_h%rrtm_aldif(1:surf_lsm_h%ns,0:2) )
+       ALLOCATE ( surf_lsm_h%rrtm_aldir(1:surf_lsm_h%ns,0:2) )
+       ALLOCATE ( surf_lsm_h%rrtm_asdif(1:surf_lsm_h%ns,0:2) )
+       ALLOCATE ( surf_lsm_h%rrtm_asdir(1:surf_lsm_h%ns,0:2) )
+       ALLOCATE ( surf_usm_h%aldif(1:surf_usm_h%ns,0:2)      )
+       ALLOCATE ( surf_usm_h%aldir(1:surf_usm_h%ns,0:2)      )
+       ALLOCATE ( surf_usm_h%asdif(1:surf_usm_h%ns,0:2)      )
+       ALLOCATE ( surf_usm_h%asdir(1:surf_usm_h%ns,0:2)      )
+       ALLOCATE ( surf_usm_h%rrtm_aldif(1:surf_usm_h%ns,0:2) )
+       ALLOCATE ( surf_usm_h%rrtm_aldir(1:surf_usm_h%ns,0:2) )
+       ALLOCATE ( surf_usm_h%rrtm_asdif(1:surf_usm_h%ns,0:2) )
+       ALLOCATE ( surf_usm_h%rrtm_asdir(1:surf_usm_h%ns,0:2) )
+!
+!--    Allocate broadband albedo (temporary for the current radiation implementations)
+       IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )                                                   &
+          ALLOCATE( surf_lsm_h%albedo(1:surf_lsm_h%ns,0:2) )
+       IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )                                                   &
+          ALLOCATE( surf_usm_h%albedo(1:surf_usm_h%ns,0:2) )
+!
+!--    Allocate albedos for short/longwave radiation, vertical surfaces
+       DO  l = 0, 3
+          ALLOCATE ( surf_lsm_v(l)%aldif(1:surf_lsm_v(l)%ns,0:2)      )
+          ALLOCATE ( surf_lsm_v(l)%aldir(1:surf_lsm_v(l)%ns,0:2)      )
+          ALLOCATE ( surf_lsm_v(l)%asdif(1:surf_lsm_v(l)%ns,0:2)      )
+          ALLOCATE ( surf_lsm_v(l)%asdir(1:surf_lsm_v(l)%ns,0:2)      )
+          ALLOCATE ( surf_lsm_v(l)%rrtm_aldif(1:surf_lsm_v(l)%ns,0:2) )
+          ALLOCATE ( surf_lsm_v(l)%rrtm_aldir(1:surf_lsm_v(l)%ns,0:2) )
+          ALLOCATE ( surf_lsm_v(l)%rrtm_asdif(1:surf_lsm_v(l)%ns,0:2) )
+          ALLOCATE ( surf_lsm_v(l)%rrtm_asdir(1:surf_lsm_v(l)%ns,0:2) )
+          ALLOCATE ( surf_usm_v(l)%aldif(1:surf_usm_v(l)%ns,0:2)      )
+          ALLOCATE ( surf_usm_v(l)%aldir(1:surf_usm_v(l)%ns,0:2)      )
+          ALLOCATE ( surf_usm_v(l)%asdif(1:surf_usm_v(l)%ns,0:2)      )
+          ALLOCATE ( surf_usm_v(l)%asdir(1:surf_usm_v(l)%ns,0:2)      )
+          ALLOCATE ( surf_usm_v(l)%rrtm_aldif(1:surf_usm_v(l)%ns,0:2) )
+          ALLOCATE ( surf_usm_v(l)%rrtm_aldir(1:surf_usm_v(l)%ns,0:2) )
+          ALLOCATE ( surf_usm_v(l)%rrtm_asdif(1:surf_usm_v(l)%ns,0:2) )
+          ALLOCATE ( surf_usm_v(l)%rrtm_asdir(1:surf_usm_v(l)%ns,0:2) )
+!
+!--       Allocate broadband albedo (temporary for the current radiation implementations)
+          IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )                                           &
+             ALLOCATE( surf_lsm_v(l)%albedo(1:surf_lsm_v(l)%ns,0:2) )
+          IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )                                           &
+             ALLOCATE( surf_usm_v(l)%albedo(1:surf_usm_v(l)%ns,0:2) )
+       ENDDO
+!
+!--    Level 1 initialization of spectral albedos via namelist paramters. Please note, this case all
+!--    surface tiles are initialized the same.
+       IF ( surf_lsm_h%ns > 0 )  THEN
+          surf_lsm_h%aldif  = albedo_lw_dif
+          surf_lsm_h%aldir  = albedo_lw_dir
+          surf_lsm_h%asdif  = albedo_sw_dif
+          surf_lsm_h%asdir  = albedo_sw_dir
+          surf_lsm_h%albedo = albedo_sw_dif
+       ENDIF
+       IF ( surf_usm_h%ns > 0 )  THEN
+          IF ( surf_usm_h%albedo_from_ascii )  THEN
+             surf_usm_h%aldif  = surf_usm_h%albedo
+             surf_usm_h%aldir  = surf_usm_h%albedo
+             surf_usm_h%asdif  = surf_usm_h%albedo
+             surf_usm_h%asdir  = surf_usm_h%albedo
+          ELSE
+             surf_usm_h%aldif  = albedo_lw_dif
+             surf_usm_h%aldir  = albedo_lw_dir
+             surf_usm_h%asdif  = albedo_sw_dif
+             surf_usm_h%asdir  = albedo_sw_dir
+             surf_usm_h%albedo = albedo_sw_dif
           ENDIF
+          IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_lw_out_change_0 )  .AND.    &
+                     surf_usm_v(l)%ns > 0  )  THEN
+             ALLOCATE( surf_usm_v(l)%rad_lw_out_change_0(1:surf_usm_v(l)%ns) )
+             surf_usm_v(l)%rad_lw_out_change_0 = 0.0_wp
+       ENDIF
+       DO  l = 0, 3
+          IF ( surf_lsm_v(l)%ns > 0 )  THEN
+             surf_lsm_v(l)%aldif  = albedo_lw_dif
+             surf_lsm_v(l)%aldir  = albedo_lw_dir
+             surf_lsm_v(l)%asdif  = albedo_sw_dif
+             surf_lsm_v(l)%asdir  = albedo_sw_dir
+             surf_lsm_v(l)%albedo = albedo_sw_dif
+          ENDIF
+          IF ( surf_usm_v(l)%ns > 0 )  THEN
+             IF ( surf_usm_v(l)%albedo_from_ascii )  THEN
+                surf_usm_v(l)%aldif  = surf_usm_v(l)%albedo
+                surf_usm_v(l)%aldir  = surf_usm_v(l)%albedo
+                surf_usm_v(l)%asdif  = surf_usm_v(l)%albedo
+                surf_usm_v(l)%asdir  = surf_usm_v(l)%albedo
+             ELSE
+                surf_usm_v(l)%aldif  = albedo_lw_dif
+                surf_usm_v(l)%aldir  = albedo_lw_dir
+                surf_usm_v(l)%asdif  = albedo_sw_dif
+                surf_usm_v(l)%asdir  = albedo_sw_dir
+             ENDIF
           ENDIF
        ENDDO
+!
+!--    Allocate surface arrays for incoming/outgoing short/longwave radiation
+       IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_sw_in )  .AND.                    &
+                  surf_lsm_h%ns > 0  )   THEN
+          ALLOCATE( surf_lsm_h%rad_sw_in(1:surf_lsm_h%ns)  )
+          ALLOCATE( surf_lsm_h%rad_sw_out(1:surf_lsm_h%ns) )
+          ALLOCATE( surf_lsm_h%rad_sw_dir(1:surf_lsm_h%ns) )
+          ALLOCATE( surf_lsm_h%rad_sw_dif(1:surf_lsm_h%ns) )
+          ALLOCATE( surf_lsm_h%rad_sw_ref(1:surf_lsm_h%ns) )
+          ALLOCATE( surf_lsm_h%rad_sw_res(1:surf_lsm_h%ns) )
+          ALLOCATE( surf_lsm_h%rad_lw_in(1:surf_lsm_h%ns)  )
+          ALLOCATE( surf_lsm_h%rad_lw_out(1:surf_lsm_h%ns) )
+          ALLOCATE( surf_lsm_h%rad_lw_dif(1:surf_lsm_h%ns) )
+          ALLOCATE( surf_lsm_h%rad_lw_ref(1:surf_lsm_h%ns) )
+          ALLOCATE( surf_lsm_h%rad_lw_res(1:surf_lsm_h%ns) )
+          surf_lsm_h%rad_sw_in  = 0.0_wp
+          surf_lsm_h%rad_sw_out = 0.0_wp
+          surf_lsm_h%rad_sw_dir = 0.0_wp
+          surf_lsm_h%rad_sw_dif = 0.0_wp
+          surf_lsm_h%rad_sw_ref = 0.0_wp
+          surf_lsm_h%rad_sw_res = 0.0_wp
+          surf_lsm_h%rad_lw_in  = 0.0_wp
+          surf_lsm_h%rad_lw_out = 0.0_wp
+          surf_lsm_h%rad_lw_dif = 0.0_wp
+          surf_lsm_h%rad_lw_ref = 0.0_wp
+          surf_lsm_h%rad_lw_res = 0.0_wp
+       ENDIF
+       IF ( .NOT. ALLOCATED ( surf_usm_h%rad_sw_in )  .AND.                    &
+                  surf_usm_h%ns > 0  )  THEN
+          ALLOCATE( surf_usm_h%rad_sw_in(1:surf_usm_h%ns)  )
+          ALLOCATE( surf_usm_h%rad_sw_out(1:surf_usm_h%ns) )
+          ALLOCATE( surf_usm_h%rad_sw_dir(1:surf_usm_h%ns) )
+          ALLOCATE( surf_usm_h%rad_sw_dif(1:surf_usm_h%ns) )
+          ALLOCATE( surf_usm_h%rad_sw_ref(1:surf_usm_h%ns) )
+          ALLOCATE( surf_usm_h%rad_sw_res(1:surf_usm_h%ns) )
+          ALLOCATE( surf_usm_h%rad_lw_in(1:surf_usm_h%ns)  )
+          ALLOCATE( surf_usm_h%rad_lw_out(1:surf_usm_h%ns) )
+          ALLOCATE( surf_usm_h%rad_lw_dif(1:surf_usm_h%ns) )
+          ALLOCATE( surf_usm_h%rad_lw_ref(1:surf_usm_h%ns) )
+          ALLOCATE( surf_usm_h%rad_lw_res(1:surf_usm_h%ns) )
+          surf_usm_h%rad_sw_in  = 0.0_wp
+          surf_usm_h%rad_sw_out = 0.0_wp
+          surf_usm_h%rad_sw_dir = 0.0_wp
+          surf_usm_h%rad_sw_dif = 0.0_wp
+          surf_usm_h%rad_sw_ref = 0.0_wp
+          surf_usm_h%rad_sw_res = 0.0_wp
+          surf_usm_h%rad_lw_in  = 0.0_wp
+          surf_usm_h%rad_lw_out = 0.0_wp
+          surf_usm_h%rad_lw_dif = 0.0_wp
+          surf_usm_h%rad_lw_ref = 0.0_wp
+          surf_usm_h%rad_lw_res = 0.0_wp
+       ENDIF
+       DO  l = 0, 3
+          IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_sw_in )  .AND.              &
+                     surf_lsm_v(l)%ns > 0  )  THEN
+             ALLOCATE( surf_lsm_v(l)%rad_sw_in(1:surf_lsm_v(l)%ns)  )
+             ALLOCATE( surf_lsm_v(l)%rad_sw_out(1:surf_lsm_v(l)%ns) )
+             ALLOCATE( surf_lsm_v(l)%rad_sw_dir(1:surf_lsm_v(l)%ns) )
+             ALLOCATE( surf_lsm_v(l)%rad_sw_dif(1:surf_lsm_v(l)%ns) )
+             ALLOCATE( surf_lsm_v(l)%rad_sw_ref(1:surf_lsm_v(l)%ns) )
+             ALLOCATE( surf_lsm_v(l)%rad_sw_res(1:surf_lsm_v(l)%ns) )
+             ALLOCATE( surf_lsm_v(l)%rad_lw_in(1:surf_lsm_v(l)%ns)  )
+             ALLOCATE( surf_lsm_v(l)%rad_lw_out(1:surf_lsm_v(l)%ns) )
+             ALLOCATE( surf_lsm_v(l)%rad_lw_dif(1:surf_lsm_v(l)%ns) )
+             ALLOCATE( surf_lsm_v(l)%rad_lw_ref(1:surf_lsm_v(l)%ns) )
+             ALLOCATE( surf_lsm_v(l)%rad_lw_res(1:surf_lsm_v(l)%ns) )
+             surf_lsm_v(l)%rad_sw_in  = 0.0_wp
+             surf_lsm_v(l)%rad_sw_out = 0.0_wp
+             surf_lsm_v(l)%rad_sw_dir = 0.0_wp
+             surf_lsm_v(l)%rad_sw_dif = 0.0_wp
+             surf_lsm_v(l)%rad_sw_ref = 0.0_wp
+             surf_lsm_v(l)%rad_sw_res = 0.0_wp
+             surf_lsm_v(l)%rad_lw_in  = 0.0_wp
+             surf_lsm_v(l)%rad_lw_out = 0.0_wp
+             surf_lsm_v(l)%rad_lw_dif = 0.0_wp
+             surf_lsm_v(l)%rad_lw_ref = 0.0_wp
+             surf_lsm_v(l)%rad_lw_res = 0.0_wp
+          ENDIF
+          IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_sw_in )  .AND.              &
+                     surf_usm_v(l)%ns > 0  )  THEN
+             ALLOCATE( surf_usm_v(l)%rad_sw_in(1:surf_usm_v(l)%ns)  )
+             ALLOCATE( surf_usm_v(l)%rad_sw_out(1:surf_usm_v(l)%ns) )
+             ALLOCATE( surf_usm_v(l)%rad_sw_dir(1:surf_usm_v(l)%ns) )
+             ALLOCATE( surf_usm_v(l)%rad_sw_dif(1:surf_usm_v(l)%ns) )
+             ALLOCATE( surf_usm_v(l)%rad_sw_ref(1:surf_usm_v(l)%ns) )
+             ALLOCATE( surf_usm_v(l)%rad_sw_res(1:surf_usm_v(l)%ns) )
+             ALLOCATE( surf_usm_v(l)%rad_lw_in(1:surf_usm_v(l)%ns)  )
+             ALLOCATE( surf_usm_v(l)%rad_lw_out(1:surf_usm_v(l)%ns) )
+             ALLOCATE( surf_usm_v(l)%rad_lw_dif(1:surf_usm_v(l)%ns) )
+             ALLOCATE( surf_usm_v(l)%rad_lw_ref(1:surf_usm_v(l)%ns) )
+             ALLOCATE( surf_usm_v(l)%rad_lw_res(1:surf_usm_v(l)%ns) )
+             surf_usm_v(l)%rad_sw_in  = 0.0_wp
+             surf_usm_v(l)%rad_sw_out = 0.0_wp
+             surf_usm_v(l)%rad_sw_dir = 0.0_wp
+             surf_usm_v(l)%rad_sw_dif = 0.0_wp
+             surf_usm_v(l)%rad_sw_ref = 0.0_wp
+             surf_usm_v(l)%rad_sw_res = 0.0_wp
+             surf_usm_v(l)%rad_lw_in  = 0.0_wp
+             surf_usm_v(l)%rad_lw_out = 0.0_wp
+             surf_usm_v(l)%rad_lw_dif = 0.0_wp
+             surf_usm_v(l)%rad_lw_ref = 0.0_wp
+             surf_usm_v(l)%rad_lw_res = 0.0_wp
+!--    Level 2 initialization of spectral albedos via albedo_type.
+!--    Please note, for natural- and urban-type surfaces, a tile approach is applied so that the
+!--    resulting albedo is calculated via the weighted average of respective surface fractions.
+       DO  m = 1, surf_lsm_h%ns
+!
+!--       Spectral albedos for vegetation/pavement/water surfaces
+          DO  ind_type = 0, 2
+             IF ( surf_lsm_h%albedo_type(m,ind_type) /= 0 )  THEN
+                surf_lsm_h%aldif(m,ind_type) = albedo_pars(1,surf_lsm_h%albedo_type(m,ind_type))
+                surf_lsm_h%asdif(m,ind_type) = albedo_pars(2,surf_lsm_h%albedo_type(m,ind_type))
+                surf_lsm_h%aldir(m,ind_type) = albedo_pars(1,surf_lsm_h%albedo_type(m,ind_type))
+                surf_lsm_h%asdir(m,ind_type) = albedo_pars(2,surf_lsm_h%albedo_type(m,ind_type))
+                surf_lsm_h%albedo(m,ind_type) = albedo_pars(0,surf_lsm_h%albedo_type(m,ind_type))
+             ENDIF
+          ENDDO
+       ENDDO
+!
+!--    For urban surface only if albedo has not already been initialized in the urban-surface model
+!--    via the ASCII file.
+       IF ( .NOT. surf_usm_h%albedo_from_ascii )  THEN
+          DO  m = 1, surf_usm_h%ns
+!
+!--          Spectral albedos for wall/green/window surfaces
+             DO  ind_type = 0, 2
+                IF ( surf_usm_h%albedo_type(m,ind_type) /= 0 )  THEN
+                   surf_usm_h%aldif(m,ind_type) = albedo_pars(1,surf_usm_h%albedo_type(m,ind_type))
+                   surf_usm_h%asdif(m,ind_type) = albedo_pars(2,surf_usm_h%albedo_type(m,ind_type))
+                   surf_usm_h%aldir(m,ind_type) = albedo_pars(1,surf_usm_h%albedo_type(m,ind_type))
+                   surf_usm_h%asdir(m,ind_type) = albedo_pars(2,surf_usm_h%albedo_type(m,ind_type))
+                   surf_usm_h%albedo(m,ind_type) = albedo_pars(0,surf_usm_h%albedo_type(m,ind_type))
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDIF
+       DO l = 0, 3
+          DO  m = 1, surf_lsm_v(l)%ns
+!
+!--          Spectral albedos for vegetation/pavement/water surfaces
+             DO  ind_type = 0, 2
+                IF ( surf_lsm_v(l)%albedo_type(m,ind_type) /= 0 )  THEN
+                   surf_lsm_v(l)%aldif(m,ind_type) =                                               &
+                         albedo_pars(1,surf_lsm_v(l)%albedo_type(m,ind_type))
+                   surf_lsm_v(l)%asdif(m,ind_type) =                                               &
+                         albedo_pars(2,surf_lsm_v(l)%albedo_type(m,ind_type))
+                   surf_lsm_v(l)%aldir(m,ind_type) =                                               &
+                         albedo_pars(1,surf_lsm_v(l)%albedo_type(m,ind_type))
+                   surf_lsm_v(l)%asdir(m,ind_type) =                                               &
+                         albedo_pars(2,surf_lsm_v(l)%albedo_type(m,ind_type))
+                   surf_lsm_v(l)%albedo(m,ind_type) =                                              &
+                         albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_type))
+                ENDIF
+             ENDDO
+          ENDDO
+!
+!--       For urban surface only if albedo has not already been initialized in the urban-surface
+!--       model via the ASCII file.
+          IF ( .NOT. surf_usm_v(l)%albedo_from_ascii )  THEN
+             DO  m = 1, surf_usm_v(l)%ns
+!
+!--             Spectral albedos for wall/green/window surfaces
+                DO  ind_type = 0, 2
+                   IF ( surf_usm_v(l)%albedo_type(m,ind_type) /= 0 )  THEN
+                      surf_usm_v(l)%aldif(m,ind_type) =                                            &
+                         albedo_pars(1,surf_usm_v(l)%albedo_type(m,ind_type))
+                      surf_usm_v(l)%asdif(m,ind_type) =                                            &
+                         albedo_pars(2,surf_usm_v(l)%albedo_type(m,ind_type))
+                      surf_usm_v(l)%aldir(m,ind_type) =                                            &
+                         albedo_pars(1,surf_usm_v(l)%albedo_type(m,ind_type))
+                      surf_usm_v(l)%asdir(m,ind_type) =                                            &
+                         albedo_pars(2,surf_usm_v(l)%albedo_type(m,ind_type))
+                      surf_usm_v(l)%albedo(m,ind_type) =                                           &
+                         albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_type))
+                   ENDIF
+                ENDDO
+             ENDDO
           ENDIF
        ENDDO
+!
+!--    Fix net radiation in case of radiation_scheme = 'constant'
+       IF ( radiation_scheme == 'constant' )  THEN
+          IF ( ALLOCATED( surf_lsm_h%rad_net ) )                               &
+             surf_lsm_h%rad_net    = net_radiation
+          IF ( ALLOCATED( surf_usm_h%rad_net ) )                               &
+             surf_usm_h%rad_net    = net_radiation
+!
+!--       Todo: weight with inclination angle
+          DO  l = 0, 3
+             IF ( ALLOCATED( surf_lsm_v(l)%rad_net ) )                         &
+                surf_lsm_v(l)%rad_net = net_radiation
+             IF ( ALLOCATED( surf_usm_v(l)%rad_net ) )                         &
+                surf_usm_v(l)%rad_net = net_radiation
+          ENDDO
+!          radiation = .FALSE.
+!
+!--    Calculate orbital constants
+       ELSE
+          decl_1 = SIN(23.45_wp * pi / 180.0_wp)
+          decl_2 = 2.0_wp * pi / 365.0_wp
+          decl_3 = decl_2 * 81.0_wp
+          lat    = latitude * pi / 180.0_wp
+          lon    = longitude * pi / 180.0_wp
+       ENDIF
+       IF ( radiation_scheme == 'clear-sky'  .OR.                              &
+            radiation_scheme == 'constant'   .OR.                              &
+            radiation_scheme == 'external' )  THEN
+!
+!--       Allocate arrays for incoming/outgoing short/longwave radiation
+          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
+             ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
+             rad_sw_in = 0.0_wp
+          ENDIF
+          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
+             ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
+             rad_sw_out = 0.0_wp
+          ENDIF
+          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
+             ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
+             rad_lw_in = 0.0_wp
+          ENDIF
+          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
+             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
+             rad_lw_out = 0.0_wp
+          ENDIF
+!
+!--       Allocate average arrays for incoming/outgoing short/longwave radiation
+          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
+             ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
+             ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
+             ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
+             ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
+          ENDIF
+!
+!--       Allocate arrays for broadband albedo, and level 1 initialization
+!--       via namelist paramter, unless not already allocated.
+          IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )  THEN
+             ALLOCATE( surf_lsm_h%albedo(1:surf_lsm_h%ns,0:2)     )
+             surf_lsm_h%albedo    = albedo
+          ENDIF
+          IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )  THEN
+             ALLOCATE( surf_usm_h%albedo(1:surf_usm_h%ns,0:2)     )
+             surf_usm_h%albedo    = albedo
+          ENDIF
+          DO  l = 0, 3
+             IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )  THEN
+                ALLOCATE( surf_lsm_v(l)%albedo(1:surf_lsm_v(l)%ns,0:2) )
+                surf_lsm_v(l)%albedo = albedo
+             ENDIF
+             IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )  THEN
+                ALLOCATE( surf_usm_v(l)%albedo(1:surf_usm_v(l)%ns,0:2) )
+                surf_usm_v(l)%albedo = albedo
+             ENDIF
+          ENDDO
+!
+!--       Level 2 initialization of broadband albedo via given albedo_type.
+!--       Only if albedo_type is non-zero. In case of urban surface and
+!--       input data is read from ASCII file, albedo_type will be zero, so that
+!--       albedo won't be overwritten.
+!--    Level 3 initialization at grid points where albedo type is zero.
+!--    In this case, spectral albedos are taken from file if available
+       IF ( albedo_pars_f%from_file )  THEN
+!
+!--       Horizontal
           DO  m = 1, surf_lsm_h%ns
+             IF ( surf_lsm_h%albedo_type(m,ind_veg_wall) /= 0 )                &
+                surf_lsm_h%albedo(m,ind_veg_wall) =                            &
+                           albedo_pars(0,surf_lsm_h%albedo_type(m,ind_veg_wall))
+             IF ( surf_lsm_h%albedo_type(m,ind_pav_green) /= 0 )               &
+                surf_lsm_h%albedo(m,ind_pav_green) =                           &
+                           albedo_pars(0,surf_lsm_h%albedo_type(m,ind_pav_green))
+             IF ( surf_lsm_h%albedo_type(m,ind_wat_win) /= 0 )                 &
+                surf_lsm_h%albedo(m,ind_wat_win) =                             &
+                           albedo_pars(0,surf_lsm_h%albedo_type(m,ind_wat_win))
+          ENDDO
+          DO  m = 1, surf_usm_h%ns
+             IF ( surf_usm_h%albedo_type(m,ind_veg_wall) /= 0 )                &
+                surf_usm_h%albedo(m,ind_veg_wall) =                            &
+                           albedo_pars(0,surf_usm_h%albedo_type(m,ind_veg_wall))
+             IF ( surf_usm_h%albedo_type(m,ind_pav_green) /= 0 )               &
+                surf_usm_h%albedo(m,ind_pav_green) =                           &
+                           albedo_pars(0,surf_usm_h%albedo_type(m,ind_pav_green))
+             IF ( surf_usm_h%albedo_type(m,ind_wat_win) /= 0 )                 &
+                surf_usm_h%albedo(m,ind_wat_win) =                             &
+                           albedo_pars(0,surf_usm_h%albedo_type(m,ind_wat_win))
+          ENDDO
+          DO  l = 0, 3
+             DO  m = 1, surf_lsm_v(l)%ns
+                IF ( surf_lsm_v(l)%albedo_type(m,ind_veg_wall) /= 0 )          &
+                   surf_lsm_v(l)%albedo(m,ind_veg_wall) =                      &
+                        albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_veg_wall))
+                IF ( surf_lsm_v(l)%albedo_type(m,ind_pav_green) /= 0 )         &
+                   surf_lsm_v(l)%albedo(m,ind_pav_green) =                     &
+                        albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_pav_green))
+                IF ( surf_lsm_v(l)%albedo_type(m,ind_wat_win) /= 0 )           &
+                   surf_lsm_v(l)%albedo(m,ind_wat_win) =                       &
+                        albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_wat_win))
+             ENDDO
+             DO  m = 1, surf_usm_v(l)%ns
+                IF ( surf_usm_v(l)%albedo_type(m,ind_veg_wall) /= 0 )          &
+                   surf_usm_v(l)%albedo(m,ind_veg_wall) =                      &
+                        albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_veg_wall))
+                IF ( surf_usm_v(l)%albedo_type(m,ind_pav_green) /= 0 )         &
+                   surf_usm_v(l)%albedo(m,ind_pav_green) =                     &
+                        albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_pav_green))
+                IF ( surf_usm_v(l)%albedo_type(m,ind_wat_win) /= 0 )           &
+                   surf_usm_v(l)%albedo(m,ind_wat_win) =                       &
+                        albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_wat_win))
+             i = surf_lsm_h%i(m)
+             j = surf_lsm_h%j(m)
+!
+!--          Spectral albedos for vegetation/pavement/water surfaces
+             DO  ind_type = 0, 2
+                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )                          &
+                   surf_lsm_h%albedo(m,ind_type) =  albedo_pars_f%pars_xy(0,j,i)
+                IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )                          &
+                   surf_lsm_h%aldir(m,ind_type) = albedo_pars_f%pars_xy(1,j,i)
+                IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )                          &
+                   surf_lsm_h%aldif(m,ind_type) = albedo_pars_f%pars_xy(1,j,i)
+                IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )                          &
+                   surf_lsm_h%asdir(m,ind_type) = albedo_pars_f%pars_xy(2,j,i)
+                IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )                          &
+                   surf_lsm_h%asdif(m,ind_type) = albedo_pars_f%pars_xy(2,j,i)
              ENDDO
           ENDDO
+!
+!--       Level 3 initialization at grid points where albedo type is zero.
+!--       This case, albedo is taken from file. In case of constant radiation
+!--       or clear sky, only broadband albedo is given.
+          IF ( albedo_pars_f%from_file )  THEN
+!
+!--          Horizontal surfaces
+             DO  m = 1, surf_lsm_h%ns
+                i = surf_lsm_h%i(m)
+                j = surf_lsm_h%j(m)
+                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
+                   surf_lsm_h%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
+                   surf_lsm_h%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
+                   surf_lsm_h%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
+                ENDIF
+             ENDDO
+!
+!--       For urban surface only if albedo has not been already initialized in the urban-surface
+!--       model via the ASCII file.
+          IF ( .NOT. surf_usm_h%albedo_from_ascii )  THEN
              DO  m = 1, surf_usm_h%ns
                 i = surf_usm_h%i(m)
                 j = surf_usm_h%j(m)
+                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
+                   surf_usm_h%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
+                   surf_usm_h%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
+                   surf_usm_h%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
+!
+!--             Broadband albedos for wall/green/window surfaces
+                DO  ind_type = 0, 2
+                   IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )                       &
+                      surf_usm_h%albedo(m,ind_type) = albedo_pars_f%pars_xy(0,j,i)
+                ENDDO
+!
+!--             Spectral albedos especially for building wall surfaces
+                IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )  THEN
+                   surf_usm_h%aldir(m,ind_veg_wall) = albedo_pars_f%pars_xy(1,j,i)
+                   surf_usm_h%aldif(m,ind_veg_wall) = albedo_pars_f%pars_xy(1,j,i)
+                ENDIF
+                IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )  THEN
+                   surf_usm_h%asdir(m,ind_veg_wall) = albedo_pars_f%pars_xy(2,j,i)
+                   surf_usm_h%asdif(m,ind_veg_wall) = albedo_pars_f%pars_xy(2,j,i)
+                ENDIF
+!
+!--             Spectral albedos especially for building green surfaces
+                IF ( albedo_pars_f%pars_xy(3,j,i) /= albedo_pars_f%fill )  THEN
+                   surf_usm_h%aldir(m,ind_pav_green) = albedo_pars_f%pars_xy(3,j,i)
+                   surf_usm_h%aldif(m,ind_pav_green) = albedo_pars_f%pars_xy(3,j,i)
+                ENDIF
+                IF ( albedo_pars_f%pars_xy(4,j,i) /= albedo_pars_f%fill )  THEN
+                   surf_usm_h%asdir(m,ind_pav_green) = albedo_pars_f%pars_xy(4,j,i)
+                   surf_usm_h%asdif(m,ind_pav_green) = albedo_pars_f%pars_xy(4,j,i)
+                ENDIF
+!
+!--             Spectral albedos especially for building window surfaces
+                IF ( albedo_pars_f%pars_xy(5,j,i) /= albedo_pars_f%fill )  THEN
+                   surf_usm_h%aldir(m,ind_wat_win) = albedo_pars_f%pars_xy(5,j,i)
+                   surf_usm_h%aldif(m,ind_wat_win) = albedo_pars_f%pars_xy(5,j,i)
+                ENDIF
+                IF ( albedo_pars_f%pars_xy(6,j,i) /= albedo_pars_f%fill )  THEN
+                   surf_usm_h%asdir(m,ind_wat_win) = albedo_pars_f%pars_xy(6,j,i)
+                   surf_usm_h%asdif(m,ind_wat_win) = albedo_pars_f%pars_xy(6,j,i)
+                ENDIF
+             ENDDO
+          ENDIF
+!
+!--       Vertical
+          DO  l = 0, 3
+             ioff = surf_lsm_v(l)%ioff
+             joff = surf_lsm_v(l)%joff
+             DO  m = 1, surf_lsm_v(l)%ns
+                i = surf_lsm_v(l)%i(m)
+                j = surf_lsm_v(l)%j(m)
+!
+!--             Spectral albedos for vegetation/pavement/water surfaces
+                DO  ind_type = 0, 2
+                   IF ( albedo_pars_f%pars_xy(0,j+joff,i+ioff) /= albedo_pars_f%fill )             &
+                      surf_lsm_v(l)%albedo(m,ind_type) = albedo_pars_f%pars_xy(0,j+joff,i+ioff)
+                   IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /= albedo_pars_f%fill )             &
+                      surf_lsm_v(l)%aldir(m,ind_type) = albedo_pars_f%pars_xy(1,j+joff,i+ioff)
+                   IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /= albedo_pars_f%fill )             &
+                      surf_lsm_v(l)%aldif(m,ind_type) = albedo_pars_f%pars_xy(1,j+joff,i+ioff)
+                   IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /= albedo_pars_f%fill )             &
+                      surf_lsm_v(l)%asdir(m,ind_type) = albedo_pars_f%pars_xy(2,j+joff,i+ioff)
+                   IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /= albedo_pars_f%fill )             &
+                      surf_lsm_v(l)%asdif(m,ind_type) = albedo_pars_f%pars_xy(2,j+joff,i+ioff)
+                ENDDO
+             ENDDO
+!
+!--          For urban surface only if albedo has not already been initialized in the urban-surface
+!--          model via the ASCII file.
+             IF ( .NOT. surf_usm_v(l)%albedo_from_ascii )  THEN
+                ioff = surf_usm_v(l)%ioff
+                joff = surf_usm_v(l)%joff
+                DO  m = 1, surf_usm_v(l)%ns
+                   i = surf_usm_v(l)%i(m)
+                   j = surf_usm_v(l)%j(m)
+!
+!--                Broadband albedos for wall/green/window surfaces
+                   DO  ind_type = 0, 2
+                      IF ( albedo_pars_f%pars_xy(0,j+joff,i+ioff) /= albedo_pars_f%fill )          &
+                         surf_usm_v(l)%albedo(m,ind_type) = albedo_pars_f%pars_xy(0,j+joff,i+ioff)
+                   ENDDO
+!
+!--                Spectral albedos especially for building wall surfaces
+                   IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /= albedo_pars_f%fill )  THEN
+                      surf_usm_v(l)%aldir(m,ind_veg_wall) = albedo_pars_f%pars_xy(1,j+joff,i+ioff)
+                      surf_usm_v(l)%aldif(m,ind_veg_wall) = albedo_pars_f%pars_xy(1,j+joff,i+ioff)
+                   ENDIF
+                   IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /= albedo_pars_f%fill )  THEN
+                      surf_usm_v(l)%asdir(m,ind_veg_wall) = albedo_pars_f%pars_xy(2,j+joff,i+ioff)
+                      surf_usm_v(l)%asdif(m,ind_veg_wall) = albedo_pars_f%pars_xy(2,j+joff,i+ioff)
+                   ENDIF
+!
+!--                Spectral albedos especially for building green surfaces
+                   IF ( albedo_pars_f%pars_xy(3,j+joff,i+ioff) /= albedo_pars_f%fill )  THEN
+                      surf_usm_v(l)%aldir(m,ind_pav_green) = albedo_pars_f%pars_xy(3,j+joff,i+ioff)
+                      surf_usm_v(l)%aldif(m,ind_pav_green) = albedo_pars_f%pars_xy(3,j+joff,i+ioff)
+                   ENDIF
+                   IF ( albedo_pars_f%pars_xy(4,j+joff,i+ioff) /= albedo_pars_f%fill )  THEN
+                      surf_usm_v(l)%asdir(m,ind_pav_green) = albedo_pars_f%pars_xy(4,j+joff,i+ioff)
+                      surf_usm_v(l)%asdif(m,ind_pav_green) = albedo_pars_f%pars_xy(4,j+joff,i+ioff)
+                   ENDIF
+!
+!--                Spectral albedos especially for building window surfaces
+                   IF ( albedo_pars_f%pars_xy(5,j+joff,i+ioff) /= albedo_pars_f%fill )  THEN
+                      surf_usm_v(l)%aldir(m,ind_wat_win) = albedo_pars_f%pars_xy(5,j+joff,i+ioff)
+                      surf_usm_v(l)%aldif(m,ind_wat_win) = albedo_pars_f%pars_xy(5,j+joff,i+ioff)
+                   ENDIF
+                   IF ( albedo_pars_f%pars_xy(6,j+joff,i+ioff) /= albedo_pars_f%fill )  THEN
+                      surf_usm_v(l)%asdir(m,ind_wat_win) = albedo_pars_f%pars_xy(6,j+joff,i+ioff)
+                      surf_usm_v(l)%asdif(m,ind_wat_win) = albedo_pars_f%pars_xy(6,j+joff,i+ioff)
+                   ENDIF
+                ENDDO
+             ENDIF
+          ENDDO
+       ENDIF
+!
+!--    Read explicit albedo values from building surface pars. If present they override all less
+!--    specific albedo values and force an albedo_type to zero in order to take effect.
+       IF ( building_surface_pars_f%from_file )  THEN
+          DO  m = 1, surf_usm_h%ns
+             i = surf_usm_h%i(m)
+             j = surf_usm_h%j(m)
+             k = surf_usm_h%k(m)
+!
+!--          Iterate over surfaces in column, check height and orientation
+             DO  is = building_surface_pars_f%index_ji(1,j,i),                                     &
+                      building_surface_pars_f%index_ji(2,j,i)
+                IF ( building_surface_pars_f%coords(4,is) == -surf_usm_h%koff  .AND.               &
+                     building_surface_pars_f%coords(1,is) == k )  THEN
+                   IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=                       &
+                        building_surface_pars_f%fill )  THEN
+                      surf_usm_h%albedo(m,ind_veg_wall) =                                          &
+                               building_surface_pars_f%pars(ind_s_alb_b_wall,is)
+                      surf_usm_h%albedo_type(m,ind_veg_wall) = 0
+                   ENDIF
+                   IF ( building_surface_pars_f%pars(ind_s_alb_l_wall,is) /=                       &
+                        building_surface_pars_f%fill )  THEN
+                      surf_usm_h%aldir(m,ind_veg_wall) =                                           &
+                               building_surface_pars_f%pars(ind_s_alb_l_wall,is)
+                      surf_usm_h%aldif(m,ind_veg_wall) =                                           &
+                               building_surface_pars_f%pars(ind_s_alb_l_wall,is)
+                      surf_usm_h%albedo_type(m,ind_veg_wall) = 0
+                   ENDIF
+                   IF ( building_surface_pars_f%pars(ind_s_alb_s_wall,is) /=                       &
+                        building_surface_pars_f%fill )  THEN
+                      surf_usm_h%asdir(m,ind_veg_wall) =                                           &
+                               building_surface_pars_f%pars(ind_s_alb_s_wall,is)
+                      surf_usm_h%asdif(m,ind_veg_wall) =                                           &
+                               building_surface_pars_f%pars(ind_s_alb_s_wall,is)
+                      surf_usm_h%albedo_type(m,ind_veg_wall) = 0
+                   ENDIF
+                   IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=                        &
+                        building_surface_pars_f%fill )  THEN
+                      surf_usm_h%albedo(m,ind_wat_win) =                                           &
+                               building_surface_pars_f%pars(ind_s_alb_b_win,is)
+                      surf_usm_h%albedo_type(m,ind_wat_win) = 0
+                   ENDIF
+                   IF ( building_surface_pars_f%pars(ind_s_alb_l_win,is) /=                        &
+                        building_surface_pars_f%fill )  THEN
+                      surf_usm_h%aldir(m,ind_wat_win) =                                            &
+                               building_surface_pars_f%pars(ind_s_alb_l_win,is)
+                      surf_usm_h%aldif(m,ind_wat_win) =                                            &
+                               building_surface_pars_f%pars(ind_s_alb_l_win,is)
+                      surf_usm_h%albedo_type(m,ind_wat_win) = 0
+                   ENDIF
+                   IF ( building_surface_pars_f%pars(ind_s_alb_s_win,is) /=                        &
+                        building_surface_pars_f%fill )  THEN
+                      surf_usm_h%asdir(m,ind_wat_win) =                                            &
+                               building_surface_pars_f%pars(ind_s_alb_s_win,is)
+                      surf_usm_h%asdif(m,ind_wat_win) =                                            &
+                               building_surface_pars_f%pars(ind_s_alb_s_win,is)
+                      surf_usm_h%albedo_type(m,ind_wat_win) = 0
+                   ENDIF
+                   IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=                      &
+                        building_surface_pars_f%fill )  THEN
+                      surf_usm_h%albedo(m,ind_pav_green) =                                         &
+                               building_surface_pars_f%pars(ind_s_alb_b_green,is)
+                      surf_usm_h%albedo_type(m,ind_pav_green) = 0
+                   ENDIF
+                   IF ( building_surface_pars_f%pars(ind_s_alb_l_green,is) /=                      &
+                        building_surface_pars_f%fill )  THEN
+                      surf_usm_h%aldir(m,ind_pav_green) =                                          &
+                               building_surface_pars_f%pars(ind_s_alb_l_green,is)
+                      surf_usm_h%aldif(m,ind_pav_green) =                                          &
+                               building_surface_pars_f%pars(ind_s_alb_l_green,is)
+                      surf_usm_h%albedo_type(m,ind_pav_green) = 0
+                   ENDIF
+                   IF ( building_surface_pars_f%pars(ind_s_alb_s_green,is) /=                      &
+                        building_surface_pars_f%fill )  THEN
+                      surf_usm_h%asdir(m,ind_pav_green) =                                          &
+                               building_surface_pars_f%pars(ind_s_alb_s_green,is)
+                      surf_usm_h%asdif(m,ind_pav_green) =                                          &
+                               building_surface_pars_f%pars(ind_s_alb_s_green,is)
+                      surf_usm_h%albedo_type(m,ind_pav_green) = 0
+                   ENDIF
+                   EXIT ! Surface was found and processed
                 ENDIF
              ENDDO
+!
+!--          Vertical surfaces
+             DO  l = 0, 3
+                ioff = surf_lsm_v(l)%ioff
+                joff = surf_lsm_v(l)%joff
+                DO  m = 1, surf_lsm_v(l)%ns
+                   i = surf_lsm_v(l)%i(m) + ioff
+                   j = surf_lsm_v(l)%j(m) + joff
+                   IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
+                      surf_lsm_v(l)%albedo(m,ind_veg_wall) = albedo_pars_f%pars_xy(0,j,i)
+                      surf_lsm_v(l)%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
+                      surf_lsm_v(l)%albedo(m,ind_wat_win) = albedo_pars_f%pars_xy(0,j,i)
+                   ENDIF
+                ENDDO
+                ioff = surf_usm_v(l)%ioff
+                joff = surf_usm_v(l)%joff
+                DO  m = 1, surf_usm_v(l)%ns
+                   i = surf_usm_v(l)%i(m) + ioff
+                   j = surf_usm_v(l)%j(m) + joff
+                   IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
+                      surf_usm_v(l)%albedo(m,ind_veg_wall) = albedo_pars_f%pars_xy(0,j,i)
+                      surf_usm_v(l)%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
+                      surf_usm_v(l)%albedo(m,ind_wat_win) = albedo_pars_f%pars_xy(0,j,i)
+          ENDDO
+          DO  l = 0, 3
+             DO  m = 1, surf_usm_v(l)%ns
+                i = surf_usm_v(l)%i(m)
+                j = surf_usm_v(l)%j(m)
+                k = surf_usm_v(l)%k(m)
+!
+!--             Iterate over surfaces in column, check height and orientation
+                DO  is = building_surface_pars_f%index_ji(1,j,i),                                  &
+                         building_surface_pars_f%index_ji(2,j,i)
+                   IF ( building_surface_pars_f%coords(5,is) == -surf_usm_v(l)%joff  .AND.         &
+                        building_surface_pars_f%coords(6,is) == -surf_usm_v(l)%ioff  .AND.         &
+                        building_surface_pars_f%coords(1,is) == k )  THEN
+                      IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=                    &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_v(l)%albedo(m,ind_veg_wall) =                                    &
+                                  building_surface_pars_f%pars(ind_s_alb_b_wall,is)
+                         surf_usm_v(l)%albedo_type(m,ind_veg_wall) = 0
+                      ENDIF
+                      IF ( building_surface_pars_f%pars(ind_s_alb_l_wall,is) /=                    &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_v(l)%aldir(m,ind_veg_wall) =                                     &
+                                  building_surface_pars_f%pars(ind_s_alb_l_wall,is)
+                         surf_usm_v(l)%aldif(m,ind_veg_wall) =                                     &
+                                  building_surface_pars_f%pars(ind_s_alb_l_wall,is)
+                         surf_usm_v(l)%albedo_type(m,ind_veg_wall) = 0
+                      ENDIF
+                      IF ( building_surface_pars_f%pars(ind_s_alb_s_wall,is) /=                    &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_v(l)%asdir(m,ind_veg_wall) =                                     &
+                                  building_surface_pars_f%pars(ind_s_alb_s_wall,is)
+                         surf_usm_v(l)%asdif(m,ind_veg_wall) =                                     &
+                                  building_surface_pars_f%pars(ind_s_alb_s_wall,is)
+                         surf_usm_v(l)%albedo_type(m,ind_veg_wall) = 0
+                      ENDIF
+                      IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=                     &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_v(l)%albedo(m,ind_wat_win) =                                     &
+                                  building_surface_pars_f%pars(ind_s_alb_b_win,is)
+                         surf_usm_v(l)%albedo_type(m,ind_wat_win) = 0
+                      ENDIF
+                      IF ( building_surface_pars_f%pars(ind_s_alb_l_win,is) /=                     &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_v(l)%aldir(m,ind_wat_win) =                                      &
+                                  building_surface_pars_f%pars(ind_s_alb_l_win,is)
+                         surf_usm_v(l)%aldif(m,ind_wat_win) =                                      &
+                                  building_surface_pars_f%pars(ind_s_alb_l_win,is)
+                         surf_usm_v(l)%albedo_type(m,ind_wat_win) = 0
+                      ENDIF
+                      IF ( building_surface_pars_f%pars(ind_s_alb_s_win,is) /=                     &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_v(l)%asdir(m,ind_wat_win) =                                      &
+                                  building_surface_pars_f%pars(ind_s_alb_s_win,is)
+                         surf_usm_v(l)%asdif(m,ind_wat_win) =                                      &
+                                  building_surface_pars_f%pars(ind_s_alb_s_win,is)
+                         surf_usm_v(l)%albedo_type(m,ind_wat_win) = 0
+                      ENDIF
+                      IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=                   &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_v(l)%albedo(m,ind_pav_green) =                                   &
+                                  building_surface_pars_f%pars(ind_s_alb_b_green,is)
+                         surf_usm_v(l)%albedo_type(m,ind_pav_green) = 0
+                      ENDIF
+                      IF ( building_surface_pars_f%pars(ind_s_alb_l_green,is) /=                   &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_v(l)%aldir(m,ind_pav_green) =                                    &
+                                  building_surface_pars_f%pars(ind_s_alb_l_green,is)
+                         surf_usm_v(l)%aldif(m,ind_pav_green) =                                    &
+                                  building_surface_pars_f%pars(ind_s_alb_l_green,is)
+                         surf_usm_v(l)%albedo_type(m,ind_pav_green) = 0
+                      ENDIF
+                      IF ( building_surface_pars_f%pars(ind_s_alb_s_green,is) /=                   &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_v(l)%asdir(m,ind_pav_green) =                                    &
+                                  building_surface_pars_f%pars(ind_s_alb_s_green,is)
+                         surf_usm_v(l)%asdif(m,ind_pav_green) =                                    &
+                                  building_surface_pars_f%pars(ind_s_alb_s_green,is)
+                         surf_usm_v(l)%albedo_type(m,ind_pav_green) = 0
+                      ENDIF
+                      EXIT ! Surface was found and processed
                    ENDIF
                 ENDDO
              ENDDO
+          ENDIF
+!
+!--       Read explicit albedo values from building surface pars. If present,
+!--       they override all less specific albedo values and force a albedo_type
+!--       to zero in order to take effect.
+          IF ( building_surface_pars_f%from_file )  THEN
+             DO  m = 1, surf_usm_h%ns
+                i = surf_usm_h%i(m)
+                j = surf_usm_h%j(m)
+                k = surf_usm_h%k(m)
+!
+!--             Iterate over surfaces in column, check height and orientation
+                DO  is = building_surface_pars_f%index_ji(1,j,i), &
+                         building_surface_pars_f%index_ji(2,j,i)
+                   IF ( building_surface_pars_f%coords(4,is) == -surf_usm_h%koff .AND. &
+                        building_surface_pars_f%coords(1,is) == k )  THEN
+                      IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=      &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_h%albedo(m,ind_veg_wall) =                         &
+                                  building_surface_pars_f%pars(ind_s_alb_b_wall,is)
+                         surf_usm_h%albedo_type(m,ind_veg_wall) = 0
+                      ENDIF
+                      IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=       &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_h%albedo(m,ind_wat_win) =                          &
+                                  building_surface_pars_f%pars(ind_s_alb_b_win,is)
+                         surf_usm_h%albedo_type(m,ind_wat_win) = 0
+                      ENDIF
+                      IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=     &
+                           building_surface_pars_f%fill )  THEN
+                         surf_usm_h%albedo(m,ind_pav_green) =                        &
+                                  building_surface_pars_f%pars(ind_s_alb_b_green,is)
+                         surf_usm_h%albedo_type(m,ind_pav_green) = 0
+                      ENDIF
+                      EXIT ! surface was found and processed
+                   ENDIF
+                ENDDO
+             ENDDO
+             DO  l = 0, 3
+                DO  m = 1, surf_usm_v(l)%ns
+                   i = surf_usm_v(l)%i(m)
+                   j = surf_usm_v(l)%j(m)
+                   k = surf_usm_v(l)%k(m)
+!
+!--                Iterate over surfaces in column, check height and orientation
+                   DO  is = building_surface_pars_f%index_ji(1,j,i), &
+                            building_surface_pars_f%index_ji(2,j,i)
+                      IF ( building_surface_pars_f%coords(5,is) == -surf_usm_v(l)%joff .AND. &
+                           building_surface_pars_f%coords(6,is) == -surf_usm_v(l)%ioff .AND. &
+                           building_surface_pars_f%coords(1,is) == k )  THEN
+                         IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=      &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_v(l)%albedo(m,ind_veg_wall) =                      &
+                                     building_surface_pars_f%pars(ind_s_alb_b_wall,is)
+                            surf_usm_v(l)%albedo_type(m,ind_veg_wall) = 0
+                         ENDIF
+                         IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=       &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_v(l)%albedo(m,ind_wat_win) =                       &
+                                     building_surface_pars_f%pars(ind_s_alb_b_win,is)
+                            surf_usm_v(l)%albedo_type(m,ind_wat_win) = 0
+                         ENDIF
+                         IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=     &
+                              building_surface_pars_f%fill )  THEN
+                            surf_usm_v(l)%albedo(m,ind_pav_green) =                     &
+                                     building_surface_pars_f%pars(ind_s_alb_b_green,is)
+                            surf_usm_v(l)%albedo_type(m,ind_pav_green) = 0
+                         ENDIF
+                         EXIT ! surface was found and processed
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+!
+!--    Initialization actions for RRTMG
+       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
+#if defined ( __rrtmg )
+!
+!--       Allocate albedos for short/longwave radiation, horizontal surfaces
+!--       for wall/green/window (USM) or vegetation/pavement/water surfaces
+!--       (LSM).
+          ALLOCATE ( surf_lsm_h%aldif(1:surf_lsm_h%ns,0:2)       )
+          ALLOCATE ( surf_lsm_h%aldir(1:surf_lsm_h%ns,0:2)       )
+          ALLOCATE ( surf_lsm_h%asdif(1:surf_lsm_h%ns,0:2)       )
+          ALLOCATE ( surf_lsm_h%asdir(1:surf_lsm_h%ns,0:2)       )
+          ALLOCATE ( surf_lsm_h%rrtm_aldif(1:surf_lsm_h%ns,0:2)  )
+          ALLOCATE ( surf_lsm_h%rrtm_aldir(1:surf_lsm_h%ns,0:2)  )
+          ALLOCATE ( surf_lsm_h%rrtm_asdif(1:surf_lsm_h%ns,0:2)  )
+          ALLOCATE ( surf_lsm_h%rrtm_asdir(1:surf_lsm_h%ns,0:2)  )
+          ALLOCATE ( surf_usm_h%aldif(1:surf_usm_h%ns,0:2)       )
+          ALLOCATE ( surf_usm_h%aldir(1:surf_usm_h%ns,0:2)       )
+          ALLOCATE ( surf_usm_h%asdif(1:surf_usm_h%ns,0:2)       )
+          ALLOCATE ( surf_usm_h%asdir(1:surf_usm_h%ns,0:2)       )
+          ALLOCATE ( surf_usm_h%rrtm_aldif(1:surf_usm_h%ns,0:2)  )
+          ALLOCATE ( surf_usm_h%rrtm_aldir(1:surf_usm_h%ns,0:2)  )
+          ALLOCATE ( surf_usm_h%rrtm_asdif(1:surf_usm_h%ns,0:2)  )
+          ALLOCATE ( surf_usm_h%rrtm_asdir(1:surf_usm_h%ns,0:2)  )
+!
+!--       Allocate broadband albedo (temporary for the current radiation
+!--       implementations)
+          IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )                            &
+             ALLOCATE( surf_lsm_h%albedo(1:surf_lsm_h%ns,0:2)     )
+          IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )                            &
+             ALLOCATE( surf_usm_h%albedo(1:surf_usm_h%ns,0:2)     )
+!
+!--       Allocate albedos for short/longwave radiation, vertical surfaces
+          ENDDO
+       ENDIF
+!
+!--    Calculate initial values of current (cosine of) zenith angle and whether the sun is up
+       CALL get_date_time( time_since_reference_point, day_of_year = day_of_year,                  &
+                           second_of_day = second_of_day )
+       CALL calc_zenith( day_of_year, second_of_day )
+!
+!--    Calculate initial surface albedo for different surfaces
+       IF ( .NOT. constant_albedo )  THEN
+#if defined( __netcdf )
+!
+!--       Horizontally aligned natural and urban surfaces
+          CALL calc_albedo( surf_lsm_h )
+          CALL calc_albedo( surf_usm_h )
+!
+!--       Vertically aligned natural and urban surfaces
           DO  l = 0, 3
+             ALLOCATE ( surf_lsm_v(l)%aldif(1:surf_lsm_v(l)%ns,0:2)      )
+             ALLOCATE ( surf_lsm_v(l)%aldir(1:surf_lsm_v(l)%ns,0:2)      )
+             ALLOCATE ( surf_lsm_v(l)%asdif(1:surf_lsm_v(l)%ns,0:2)      )
+             ALLOCATE ( surf_lsm_v(l)%asdir(1:surf_lsm_v(l)%ns,0:2)      )
+             ALLOCATE ( surf_lsm_v(l)%rrtm_aldif(1:surf_lsm_v(l)%ns,0:2) )
+             ALLOCATE ( surf_lsm_v(l)%rrtm_aldir(1:surf_lsm_v(l)%ns,0:2) )
+             ALLOCATE ( surf_lsm_v(l)%rrtm_asdif(1:surf_lsm_v(l)%ns,0:2) )
+             ALLOCATE ( surf_lsm_v(l)%rrtm_asdir(1:surf_lsm_v(l)%ns,0:2) )
+             ALLOCATE ( surf_usm_v(l)%aldif(1:surf_usm_v(l)%ns,0:2)      )
+             ALLOCATE ( surf_usm_v(l)%aldir(1:surf_usm_v(l)%ns,0:2)      )
+             ALLOCATE ( surf_usm_v(l)%asdif(1:surf_usm_v(l)%ns,0:2)      )
+             ALLOCATE ( surf_usm_v(l)%asdir(1:surf_usm_v(l)%ns,0:2)      )
+             ALLOCATE ( surf_usm_v(l)%rrtm_aldif(1:surf_usm_v(l)%ns,0:2) )
+             ALLOCATE ( surf_usm_v(l)%rrtm_aldir(1:surf_usm_v(l)%ns,0:2) )
+             ALLOCATE ( surf_usm_v(l)%rrtm_asdif(1:surf_usm_v(l)%ns,0:2) )
+             ALLOCATE ( surf_usm_v(l)%rrtm_asdir(1:surf_usm_v(l)%ns,0:2) )
+!
+!--          Allocate broadband albedo (temporary for the current radiation
+!--          implementations)
+             IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )                    &
+                ALLOCATE( surf_lsm_v(l)%albedo(1:surf_lsm_v(l)%ns,0:2) )
+             IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )                    &
+                ALLOCATE( surf_usm_v(l)%albedo(1:surf_usm_v(l)%ns,0:2) )
+             CALL calc_albedo( surf_lsm_v(l) )
+             CALL calc_albedo( surf_usm_v(l) )
           ENDDO
+!
+!--       Level 1 initialization of spectral albedos via namelist
+!--       paramters. Please note, this case all surface tiles are initialized
+!--       the same.
+#endif
+       ELSE
+!
+!--       Initialize sun-inclination independent spectral albedos
+!--       Horizontal surfaces
           IF ( surf_lsm_h%ns > 0 )  THEN
+             surf_lsm_h%aldif  = albedo_lw_dif
+             surf_lsm_h%aldir  = albedo_lw_dir
+             surf_lsm_h%asdif  = albedo_sw_dif
+             surf_lsm_h%asdir  = albedo_sw_dir
+    &nbs