source: palm/trunk/SOURCE/radiation_model_mod.f90 @ 4783

!> @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
! Copyright 2015-2020 Czech Technical University in Prague
! Copyright 1997-2020 Leibniz Universitaet Hannover
!------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
!
!
! Former revisions:
! -----------------
! $Id: radiation_model_mod.f90 4783 2020-11-13 13:58:45Z raasch $
! Albedo parameter list extended - building facades added
!
! 4780 2020-11-10 11:17:10Z suehring
! Enable 3D data output also with 64-bit precision
!
! 4741 2020-10-14 14:32:50Z suehring
! Add option to force calculation of horizontal mean profiles independent on data output
!
! 4717 2020-09-30 22:27:40Z pavelkrc
! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
! directives (J. Resler)
!
! 4713 2020-09-29 12:02:05Z pavelkrc
! Correct OpenMP parallelization including cycles with cumulative variables (J. Resler)
!
! 4708 2020-09-28 17:42:58Z suehring
! - Bugfix, correct mapping of RRTMG heating rates, as well as incoming/outgoing
!   radiation onto the topography-following grid
! - add fill values to the output
!
! 4694 2020-09-23 15:09:19Z pavelkrc
! Bugfix for tracing of maximum radiative fluxes
!
! 4686 2020-09-18 11:59:57Z pavelkrc
! Silence compiler warning about unused variable in serial build
!
! 4684 2020-09-18 09:50:28Z pavelkrc
! Bugfix to previous for serial build
!
! 4683 2020-09-18 09:23:09Z pavelkrc
! Add option to limit the size of MPI_Alltoall calls
!
! 4679 2020-09-14 13:53:52Z pavelkrc
! Enable warning about unrealistically large radiative fluxes by default
!
! 4678 2020-09-14 10:23:14Z pavelkrc
! Fix zdirs indices for full3d in raytrace_2d
!
! 4671 2020-09-09 20:27:58Z pavelkrc
! Radiative transfer model RTM version 4.1
! - Implementation of downward facing USM and LSM surfaces
! - Removal of deprecated CSV inputs
! - Bugfixes
! Authors: J. Resler, P. Krc  (Institute of Computer Science, Prague)
!
! 4668 2020-09-09 13:00:16Z pavelkrc
! Improve debug messages during timestepping
!
! 4664 2020-09-02 15:36:10Z gronemeier
! Add additional debug messages around MPI-AlltoAllv calls
!
! 4662 2020-09-02 13:40:38Z pavelkrc
! Bugfix of LAD in CSF generation in legacy raytracing
!
! 4661 2020-09-01 16:00:22Z pavelkrc
! Bugfix for rad_angular_discretization = .FALSE.
!
! 4657 2020-08-28 16:53:45Z pavelkrc
! Bugfix (silence compiler warning)
!
! 4656 2020-08-28 14:51:30Z pavelkrc
! Bugfix for serial build
!
! 4653 2020-08-27 08:54:43Z pavelkrc
! RTM version 4.0:
! - Support full-3D geometry:
!   * Downward faces can only be default type (temporary PALM limitation)
!   * rad_angular_discretization = .FALSE. is not supported for full-3D cases
! - Order faces in surf(l) by column
! - Remove separate direction indices for urban and land faces
! - Bugfix: skyvft output
! Authors: P. Krc, J. Resler (Institute of Computer Science, Prague)
!
! 4644 2020-08-22 15:05:46Z pavelkrc
! Bugfix: remove invalid reordering of vertical surfaces in RTM
!
! 4643 2020-08-17 09:19:42Z pavelkrc
! - Bugfix: typo in index of surf%* arrays
!
! 4634 2020-08-05 15:42:28Z pavelkrc
! - Bugfix: rtm_svf_, rtm_dif_ outputs
! - Bugfix: correct average transparency for MRT factors
!
! 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
!
! 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
!
! 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
!
! 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
!
! 4517 2020-05-03 14:29:30Z raasch
! added restart with MPI-IO for reading local arrays
!
! 4495 2020-04-13 20:11:20Z raasch
! 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
!
! 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.
! - 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
!
! 4400 2020-02-10 20:32:41Z suehring
! Initialize radiation arrays with zero
!
! 4392 2020-01-31 16:14:57Z pavelkrc
! - Add debug tracing of large radiative fluxes (option trace_fluxes_above)
! - Print exact counts of SVF and CSF if debut_output is enabled
! - Update descriptions of RTM 3.0 and related comments
!
! 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
!
! 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
!
! 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
!
! 4286 2019-10-30 16:01:14Z resler
! - Fix wrong treating of time_rad during interpolation in radiation_model_mod
! - Fix wrong checks of time_rad from dynamic driver in radiation_model_mod
! - Add new directional model of human body for MRT: ellipsoid
!
! 4271 2019-10-23 10:46:41Z maronga
! Bugfix: missing parentheses in calculation of snow albedo
!
! 4245 2019-09-30 08:40:37Z pavelkrc
! Initialize explicit per-surface albedos from building_surface_pars
!
! 4238 2019-09-25 16:06:01Z suehring
! Modify check in order to avoid equality comparisons of floating points
!
! 4227 2019-09-10 18:04:34Z gronemeier
! implement new palm_date_time_mod
!
! 4226 2019-09-10 17:03:24Z suehring
! - Netcdf input routine for dimension length renamed
! - Define time variable for external radiation input relative to time_utc_init
!
! 4210 2019-09-02 13:07:09Z suehring
! - Revise steering of splitting diffuse and direct radiation
! - Bugfixes in checks
! - 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)
!
! 4198 2019-08-29 15:17:48Z gronemeier
! Prohibit execution of radiation model if rotation_angle is not zero
!
! 4197 2019-08-29 14:33:32Z suehring
! Revise steering of surface albedo initialization when albedo_pars is provided
!
! 4190 2019-08-27 15:42:37Z suehring
! Implement external radiation forcing also for level-of-detail = 2
! (horizontally 2D radiation)
!
! 4188 2019-08-26 14:15:47Z suehring
! Minor adjustment in error message
!
! 4187 2019-08-26 12:43:15Z suehring
! - Take external radiation from root domain dynamic input if not provided for
!   each nested domain
! - Combine MPI_ALLREDUCE calls to reduce mpi overhead
!
! 4182 2019-08-22 15:20:23Z scharf
! Corrected "Former revisions" section
!
! 4179 2019-08-21 11:16:12Z suehring
! Remove debug prints
!
! 4178 2019-08-21 11:13:06Z suehring
! External radiation forcing implemented.
!
! 4168 2019-08-16 13:50:17Z suehring
! Replace function get_topography_top_index by topo_top_ind
!
! 4157 2019-08-14 09:19:12Z suehring
! Give informative message on raytracing distance only by core zero
!
! 4148 2019-08-08 11:26:00Z suehring
! Comments added
!
! 4134 2019-08-02 18:39:57Z suehring
! Bugfix in formatted write statement
!
! 4127 2019-07-30 14:47:10Z suehring
! Remove unused pch_index (merge from branch resler)
!
! 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
!
! 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
!
! 4067 2019-07-01 13:29:25Z suehring
! Bugfix, pass dummy string to MPI_INFO_SET (J. Resler)
!
! 4039 2019-06-18 10:32:41Z suehring
! Bugfix for masked data output
!
! 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
!
! 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
!
! 3987 2019-05-22 09:52:13Z kanani
! Introduce alternative switch for debug output during timestepping
!
! 3943 2019-05-02 09:50:41Z maronga
! Missing blank characteer added.
!
! 3900 2019-04-16 15:17:43Z suehring
! Fixed initialization problem
!
! 3885 2019-04-11 11:29:34Z kanani
! 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
!
! 3861 2019-04-04 06:27:41Z maronga
! Bugfix: factor of 4.0 required instead of 3.0 in calculation of rad_lw_out_change_0
!
! 3859 2019-04-03 20:30:31Z maronga
! Added some descriptions
!
! 3847 2019-04-01 14:51:44Z suehring
! Implement check for dt_radiation (must be > 0)
!
! 3846 2019-04-01 13:55:30Z suehring
! unused variable removed
!
! 3814 2019-03-26 08:40:31Z pavelkrc
! Change zenith(0:0) and others to scalar.
! Code review.
! Rename exported nzu, nzp and related variables due to name conflict
!
! 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
!
! 3769 2019-02-28 10:16:49Z moh.hefny
! 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
!
! 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).
!
! 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.
!
! 3752 2019-02-19 09:37:22Z resler
! added read/write number of MRT factors to the respective routines
!
! 3705 2019-01-29 19:56:39Z suehring
! Make variables that are sampled in virtual measurement module public
!
! 3704 2019-01-29 19:51:41Z suehring
! Some interface calls moved to module_interface + cleanup
!
! 3667 2019-01-10 14:26:24Z schwenkel
! Modified check for rrtmg input files
!
! 3655 2019-01-07 16:51:22Z knoop
! nopointer option removed
!
! 1496 2014-12-02 17:25:50Z maronga
! Initial revision
!
!
! Description:
! ------------
!> Radiation models and interfaces:
!> 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
!> Integrations of RTM with other PALM-4U 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 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.
!------------------------------------------------------------------------------!
 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,                                                 &
        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, wall_flags_total_0

    USE, INTRINSIC :: iso_c_binding

    USE kinds

    USE bulk_cloud_model_mod,                                                  &
        ONLY:  bulk_cloud_model, microphysics_morrison, na_init, nc_const, sigma_gc

#if defined ( __netcdf )
    USE 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 pegrid

#if defined ( __rrtmg )
    USE parrrsw,                                                               &
        ONLY:  naerec, nbndsw

    USE parrrtm,                                                               &
        ONLY:  nbndlw

    USE rrtmg_lw_init,                                                         &
        ONLY:  rrtmg_lw_ini

    USE rrtmg_sw_init,                                                         &
        ONLY:  rrtmg_sw_ini

    USE rrtmg_lw_rad,                                                          &
        ONLY:  rrtmg_lw

    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:  hom

    USE surface_mod,                                                           &
        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h,           &
               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'

!
!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
    CHARACTER(37), DIMENSION(0:42), 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
                                   'building wall - reflecting facade    ', & ! 34
                                   'building wall - bright facade        ', & ! 35
                                   'building wall - other materials      ', & ! 36
                                   'building window - double glazing     ', & ! 37
                                   'building window - double glazing     ', & ! 38
                                   'building window - reflecting         ', & ! 39
                                   'building roof - reflecting           ', & ! 40
                                   'building roof - bright               ', & ! 41
                                   'building roof - other materials      '  & ! 42
                                                         /)
!
!-- 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 = 2000.0_wp,  & !< NAMELIST option for debug printing of largest radiative fluxes
                                                   !< (W/m2 for surfaces, W/m3 for PC). -1=off, 0=all fluxes
                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
!
!-- Land surface albedos for solar zenith angle of 60degree after Briegleb (1992)
!-- (broadband, longwave, shortwave ):   bb,      lw,      sw,
    REAL(wp), DIMENSION(0:2,1:42), PARAMETER :: albedo_pars = RESHAPE( (/&
                                   0.06_wp, 0.06_wp, 0.06_wp,            & !  1 - ocean
                                   0.19_wp, 0.28_wp, 0.09_wp,            & !  2 - mixed farming, tall grassland
                                   0.23_wp, 0.33_wp, 0.11_wp,            & !  3 - tall/medium grassland 
                                   0.23_wp, 0.33_wp, 0.11_wp,            & !  4 - evergreen shrubland 
                                   0.25_wp, 0.34_wp, 0.14_wp,            & !  5 - short grassland/meadow/shrubland 
                                   0.14_wp, 0.22_wp, 0.06_wp,            & !  6 - evergreen needleleaf forest 
                                   0.17_wp, 0.27_wp, 0.06_wp,            & !  7 - mixed deciduous forest 
                                   0.19_wp, 0.31_wp, 0.06_wp,            & !  8 - deciduous forest 
                                   0.14_wp, 0.22_wp, 0.06_wp,            & !  9 - tropical evergreen broadleaved forest
                                   0.18_wp, 0.28_wp, 0.06_wp,            & ! 10 - medium/tall grassland/woodland 
                                   0.43_wp, 0.51_wp, 0.35_wp,            & ! 11 - desert, sandy 
                                   0.32_wp, 0.40_wp, 0.24_wp,            & ! 12 - desert, rocky 
                                   0.19_wp, 0.27_wp, 0.10_wp,            & ! 13 - tundra 
                                   0.77_wp, 0.65_wp, 0.90_wp,            & ! 14 - land ice 
                                   0.77_wp, 0.65_wp, 0.90_wp,            & ! 15 - sea ice 
                                   0.82_wp, 0.70_wp, 0.95_wp,            & ! 16 - snow 
                                   0.08_wp, 0.08_wp, 0.08_wp,            & ! 17 - bare soil 
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 18 - asphalt/concrete mix 
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 19 - asphalt (asphalt concrete) 
                                   0.30_wp, 0.30_wp, 0.30_wp,            & ! 20 - concrete (Portland concrete) 
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 21 - sett 
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 22 - paving stones
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 23 - cobblestone 
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 24 - metal
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 25 - wood
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 26 - gravel
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 27 - fine gravel 
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 28 - pebblestone 
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 29 - woodchips
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 30 - tartan (sports)
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 31 - artificial turf (sports)
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 32 - clay (sports) 
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 33 - building (dummy)
                                   0.60_wp, 0.60_wp, 0.60_wp,            & ! 34 - building wall type 1 - reflecting facade
                                   0.30_wp, 0.30_wp, 0.30_wp,            & ! 35 - building wall type 2 - bright facacde
                                   0.07_wp, 0.07_wp, 0.07_wp,            & ! 36 - building wall type 3) - other materials
                                   0.12_wp, 0.12_wp, 0.12_wp,            & ! 37 - building window type 1 - double glazing
                                   0.17_wp, 0.18_wp, 0.18_wp,            & ! 38 - building window type 2 - triple glazing
                                   0.48_wp, 0.48_wp, 0.48_wp,            & ! 39 - building window type 3 - reflecting
                                   0.60_wp, 0.60_wp, 0.60_wp,            & ! 40 - building roof type 1 - reflecting
                                   0.30_wp, 0.30_wp, 0.30_wp,            & ! 41 - building roof type 2 - bright
                                   0.07_wp, 0.07_wp, 0.07_wp             & ! 42 - building roof type 3 - other materials
                                 /), (/ 3, 42 /) )

    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)


!
!-- 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

!
!-- 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)

#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                        ::  nidx_surf = 4                      !< number of indices in surfl and surf

    INTEGER(iwp), PARAMETER                        ::  nsurf_type = 5                     !< number of surf types = surface directions

    INTEGER(iwp), PARAMETER                        ::  iup    = 0                         !< 0 - index of upward surface (ground or roof)
    INTEGER(iwp), PARAMETER                        ::  idown  = 1                         !< 1 - index of downward surface (overhanging)
    INTEGER(iwp), PARAMETER                        ::  inorth = 2                         !< 2 - index of northward facing wall
    INTEGER(iwp), PARAMETER                        ::  isouth = 3                         !< 3 - index of southward facing wall
    INTEGER(iwp), PARAMETER                        ::  ieast  = 4                         !< 4 - index of eastward facing wall
    INTEGER(iwp), PARAMETER                        ::  iwest  = 5                         !< 5 - index of westward facing wall

    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  idir = (/0, 0, 0, 0, 1,-1/)      !< surface normal direction x indices
    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  jdir = (/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/)      !< surface normal direction z indices
    REAL(wp),     DIMENSION(0:nsurf_type)          :: facearea                            !< area of single face in respective
                                                                                          !< direction (will be calc'd)


!-- indices needed for RTM netcdf output subroutines
    INTEGER(iwp), PARAMETER                        ::  nd = 6                                  !< number of directions
    CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER ::  dirname = (/ '_up   ', '_down ', '_south', '_north', '_west ', '_east ' /)
    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     ::  dirint = (/ iup, idown, isouth, inorth, iwest, ieast /) !< direction integers
    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     ::  diridx =  (/  0,     1,      1,      0,     3,     2 /) !< mapping to surf_h and surf_v
    INTEGER(iwp), DIMENSION(0:1), PARAMETER        ::  dirint_h = (/ iup, idown /)                             !< mapping to surf_h directions
    INTEGER(iwp), DIMENSION(0:3), PARAMETER        ::  dirint_v = (/ inorth, isouth, ieast, iwest /)           !< mapping to surf_v directions

!-- 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]
    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_col_start  !< start of surfaces in surfl for each x,y column (local surfaces)
    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfg_col_start  !< start of surfaces in surfl for each x,y column (all surfaces)

!-- 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)
    INTEGER(iwp)                                   ::  bufsize_alltoall = 0               !< max no. of items to send in mpi_alltoall at once (0=infinite)
    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. 4.1'         !< 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
    TYPE t_svf
        INTEGER(iwp)                               :: isurflt           !<
        INTEGER(iwp)                               :: isurfs            !<
        REAL(wp)                                   :: rsvf              !<
        REAL(wp)                                   :: rtransp           !<
    END TYPE

!-- 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)
    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
#if defined( __parallel )
    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        ::  nzterrt          !< temporary global arrays for raytracing
    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nzterrb          !< temporary global arrays for raytracing
    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  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
#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
#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


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-- 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

    INTERFACE radiation_check_data_output
       MODULE PROCEDURE radiation_check_data_output
    END INTERFACE radiation_check_data_output

    INTERFACE radiation_check_data_output_ts
       MODULE PROCEDURE radiation_check_data_output_ts
    END INTERFACE radiation_check_data_output_ts

    INTERFACE radiation_check_data_output_pr
       MODULE PROCEDURE radiation_check_data_output_pr
    END INTERFACE radiation_check_data_output_pr

    INTERFACE radiation_check_parameters
       MODULE PROCEDURE radiation_check_parameters
    END INTERFACE radiation_check_parameters

    INTERFACE radiation_clearsky
       MODULE PROCEDURE radiation_clearsky
    END INTERFACE radiation_clearsky

    INTERFACE radiation_constant
       MODULE PROCEDURE radiation_constant
    END INTERFACE radiation_constant

    INTERFACE radiation_control
       MODULE PROCEDURE radiation_control
    END INTERFACE radiation_control

    INTERFACE radiation_3d_data_averaging
       MODULE PROCEDURE radiation_3d_data_averaging
    END INTERFACE radiation_3d_data_averaging

    INTERFACE radiation_data_output_2d
       MODULE PROCEDURE radiation_data_output_2d
    END INTERFACE radiation_data_output_2d

    INTERFACE radiation_data_output_3d
       MODULE PROCEDURE radiation_data_output_3d
    END INTERFACE radiation_data_output_3d

    INTERFACE radiation_data_output_mask
       MODULE PROCEDURE radiation_data_output_mask
    END INTERFACE radiation_data_output_mask

    INTERFACE radiation_define_netcdf_grid
       MODULE PROCEDURE radiation_define_netcdf_grid
    END INTERFACE radiation_define_netcdf_grid

    INTERFACE radiation_header
       MODULE PROCEDURE radiation_header
    END INTERFACE radiation_header

    INTERFACE radiation_init
       MODULE PROCEDURE radiation_init
    END INTERFACE radiation_init

    INTERFACE radiation_parin
       MODULE PROCEDURE radiation_parin
    END INTERFACE radiation_parin

    INTERFACE radiation_rrtmg
       MODULE PROCEDURE radiation_rrtmg
    END INTERFACE radiation_rrtmg

#if defined( __rrtmg )
    INTERFACE radiation_tendency
       MODULE PROCEDURE radiation_tendency
       MODULE PROCEDURE radiation_tendency_ij
    END INTERFACE radiation_tendency
#endif

    INTERFACE radiation_rrd_local
       MODULE PROCEDURE radiation_rrd_local_ftn
       MODULE PROCEDURE radiation_rrd_local_mpi
    END INTERFACE radiation_rrd_local

    INTERFACE radiation_wrd_local
       MODULE PROCEDURE radiation_wrd_local
    END INTERFACE radiation_wrd_local

    INTERFACE radiation_interaction
       MODULE PROCEDURE radiation_interaction
    END INTERFACE radiation_interaction

    INTERFACE radiation_interaction_init
       MODULE PROCEDURE radiation_interaction_init
    END INTERFACE radiation_interaction_init

    INTERFACE radiation_presimulate_solar_pos
       MODULE PROCEDURE radiation_presimulate_solar_pos
    END INTERFACE radiation_presimulate_solar_pos

    INTERFACE radiation_calc_svf
       MODULE PROCEDURE radiation_calc_svf
    END INTERFACE radiation_calc_svf

    INTERFACE radiation_write_svf
       MODULE PROCEDURE radiation_write_svf
    END INTERFACE radiation_write_svf

    INTERFACE radiation_read_svf
       MODULE PROCEDURE radiation_read_svf
    END INTERFACE radiation_read_svf


    SAVE

    PRIVATE

!
!-- 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 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, idown, inorth, isouth, ieast, iwest,                           &
           nd, dirname, diridx, dirint,                                        &
           nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf,           &
           idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
           radiation_interactions, skyvf, skyvft, radiation_interactions_on,   &
           average_radiation, rad_sw_in_diff, rad_sw_in_dir


#if defined ( __rrtmg )
    PUBLIC radiation_tendency, rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
#endif

 CONTAINS


!------------------------------------------------------------------------------!
! 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' )
             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

!------------------------------------------------------------------------------!
! 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
       INTEGER(iwp)                 ::  id

       var = TRIM(variable)
!
!--    Identify directional variables
       ilast_word = SCAN(var, '_', back=.TRUE.)
       directional = .FALSE.
       IF ( ilast_word > 0 )  THEN
          DO id = 0, nd-1
             IF ( TRIM(var(ilast_word:)) == TRIM(dirname(id)) ) THEN
                directional = .TRUE.
                var = var(1:ilast_word-1)
                EXIT
             ENDIF
          ENDDO
       ENDIF

       IF ( directional )  THEN
          IF ( var(1:8) == 'rtm_svf_'  .OR.  var(1:8) == 'rtm_dif_' )  THEN
             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'
          ELSE
             SELECT CASE ( TRIM(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_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
          ENDIF

       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


!------------------------------------------------------------------------------!
! 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

!
!-- Next line is just to avoid compiler warning about unused variable.
    IF ( dots_max == 0 )  CONTINUE

!
!-- 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
          dots_num = dots_num + 15
       ELSE
          dots_num = dots_num + 11
       ENDIF
    ENDIF


 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


!------------------------------------------------------------------------------!
! 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
#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 )
#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 )
#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
#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
#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
#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


!------------------------------------------------------------------------------!
! 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
#if defined( __rrtmg )
       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
          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
       DO  l = 0, 1
          IF ( .NOT. ALLOCATED ( surf_lsm_h(l)%rad_net )  .AND.                      &
                     surf_lsm_h(l)%ns > 0  )   THEN
             ALLOCATE( surf_lsm_h(l)%rad_net(1:surf_lsm_h(l)%ns) )
             surf_lsm_h(l)%rad_net = 0.0_wp
          ENDIF
          IF ( .NOT. ALLOCATED ( surf_usm_h(l)%rad_net )  .AND.                      &
                     surf_usm_h(l)%ns > 0  )  THEN
             ALLOCATE( surf_usm_h(l)%rad_net(1:surf_usm_h(l)%ns) )
             surf_usm_h(l)%rad_net = 0.0_wp
          ENDIF
       ENDDO
       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
       DO  l = 0, 1
          IF ( .NOT. ALLOCATED ( surf_lsm_h(l)%rad_lw_out_change_0 )  .AND.          &
                     surf_lsm_h(l)%ns > 0  )   THEN
             ALLOCATE( surf_lsm_h(l)%rad_lw_out_change_0(1:surf_lsm_h(l)%ns) )
             surf_lsm_h(l)%rad_lw_out_change_0 = 0.0_wp
          ENDIF
          IF ( .NOT. ALLOCATED ( surf_usm_h(l)%rad_lw_out_change_0 )  .AND.          &
                     surf_usm_h(l)%ns > 0  )  THEN
             ALLOCATE( surf_usm_h(l)%rad_lw_out_change_0(1:surf_usm_h(l)%ns) )
             surf_usm_h(l)%rad_lw_out_change_0 = 0.0_wp
          ENDIF
       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
          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
       DO  l = 0, 1
          IF ( .NOT. ALLOCATED ( surf_lsm_h(l)%rad_sw_in )  .AND.                    &
                     surf_lsm_h(l)%ns > 0  )   THEN
             ALLOCATE( surf_lsm_h(l)%rad_sw_in(1:surf_lsm_h(l)%ns)  )
             ALLOCATE( surf_lsm_h(l)%rad_sw_out(1:surf_lsm_h(l)%ns) )
             ALLOCATE( surf_lsm_h(l)%rad_sw_dir(1:surf_lsm_h(l)%ns) )
             ALLOCATE( surf_lsm_h(l)%rad_sw_dif(1:surf_lsm_h(l)%ns) )
             ALLOCATE( surf_lsm_h(l)%rad_sw_ref(1:surf_lsm_h(l)%ns) )
             ALLOCATE( surf_lsm_h(l)%rad_sw_res(1:surf_lsm_h(l)%ns) )
             ALLOCATE( surf_lsm_h(l)%rad_lw_in(1:surf_lsm_h(l)%ns)  )
             ALLOCATE( surf_lsm_h(l)%rad_lw_out(1:surf_lsm_h(l)%ns) )
             ALLOCATE( surf_lsm_h(l)%rad_lw_dif(1:surf_lsm_h(l)%ns) )
             ALLOCATE( surf_lsm_h(l)%rad_lw_ref(1:surf_lsm_h(l)%ns) )
             ALLOCATE( surf_lsm_h(l)%rad_lw_res(1:surf_lsm_h(l)%ns) )
             surf_lsm_h(l)%rad_sw_in  = 0.0_wp
             surf_lsm_h(l)%rad_sw_out = 0.0_wp
             surf_lsm_h(l)%rad_sw_dir = 0.0_wp
             surf_lsm_h(l)%rad_sw_dif = 0.0_wp
             surf_lsm_h(l)%rad_sw_ref = 0.0_wp
             surf_lsm_h(l)%rad_sw_res = 0.0_wp
             surf_lsm_h(l)%rad_lw_in  = 0.0_wp
             surf_lsm_h(l)%rad_lw_out = 0.0_wp
             surf_lsm_h(l)%rad_lw_dif = 0.0_wp
             surf_lsm_h(l)%rad_lw_ref = 0.0_wp
             surf_lsm_h(l)%rad_lw_res = 0.0_wp
          ENDIF
          IF ( .NOT. ALLOCATED ( surf_usm_h(l)%rad_sw_in )  .AND.                    &
                     surf_usm_h(l)%ns > 0  )  THEN
             ALLOCATE( surf_usm_h(l)%rad_sw_in(1:surf_usm_h(l)%ns)  )
             ALLOCATE( surf_usm_h(l)%rad_sw_out(1:surf_usm_h(l)%ns) )
             ALLOCATE( surf_usm_h(l)%rad_sw_dir(1:surf_usm_h(l)%ns) )
             ALLOCATE( surf_usm_h(l)%rad_sw_dif(1:surf_usm_h(l)%ns) )
             ALLOCATE( surf_usm_h(l)%rad_sw_ref(1:surf_usm_h(l)%ns) )
             ALLOCATE( surf_usm_h(l)%rad_sw_res(1:surf_usm_h(l)%ns) )
             ALLOCATE( surf_usm_h(l)%rad_lw_in(1:surf_usm_h(l)%ns)  )
             ALLOCATE( surf_usm_h(l)%rad_lw_out(1:surf_usm_h(l)%ns) )
             ALLOCATE( surf_usm_h(l)%rad_lw_dif(1:surf_usm_h(l)%ns) )
             ALLOCATE( surf_usm_h(l)%rad_lw_ref(1:surf_usm_h(l)%ns) )
             ALLOCATE( surf_usm_h(l)%rad_lw_res(1:surf_usm_h(l)%ns) )
             surf_usm_h(l)%rad_sw_in  = 0.0_wp
             surf_usm_h(l)%rad_sw_out = 0.0_wp
             surf_usm_h(l)%rad_sw_dir = 0.0_wp
             surf_usm_h(l)%rad_sw_dif = 0.0_wp
             surf_usm_h(l)%rad_sw_ref = 0.0_wp
             surf_usm_h(l)%rad_sw_res = 0.0_wp
             surf_usm_h(l)%rad_lw_in  = 0.0_wp
             surf_usm_h(l)%rad_lw_out = 0.0_wp
             surf_usm_h(l)%rad_lw_dif = 0.0_wp
             surf_usm_h(l)%rad_lw_ref = 0.0_wp
             surf_usm_h(l)%rad_lw_res = 0.0_wp
          ENDIF
       ENDDO
       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
          DO  l = 0, 1
             IF ( ALLOCATED( surf_lsm_h(l)%rad_net ) )                               &
                surf_lsm_h(l)%rad_net    = net_radiation
             IF ( ALLOCATED( surf_usm_h(l)%rad_net ) )                               &
                surf_usm_h(l)%rad_net    = net_radiation
          ENDDO
!
!--       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.
          DO  l = 0, 1
             IF ( .NOT. ALLOCATED(surf_lsm_h(l)%albedo) )  THEN
                ALLOCATE( surf_lsm_h(l)%albedo(1:surf_lsm_h(l)%ns,0:2)     )
                surf_lsm_h(l)%albedo    = albedo
             ENDIF
             IF ( .NOT. ALLOCATED(surf_usm_h(l)%albedo) )  THEN
                ALLOCATE( surf_usm_h(l)%albedo(1:surf_usm_h(l)%ns,0:2)     )
                surf_usm_h(l)%albedo    = albedo
             ENDIF
          ENDDO
          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.
          DO  l = 0, 1
             DO  m = 1, surf_lsm_h(l)%ns
                IF ( surf_lsm_h(l)%albedo_type(m,ind_veg_wall) /= 0 )                &
                   surf_lsm_h(l)%albedo(m,ind_veg_wall) =                            &
                              albedo_pars(0,surf_lsm_h(l)%albedo_type(m,ind_veg_wall))
                IF ( surf_lsm_h(l)%albedo_type(m,ind_pav_green) /= 0 )               &
                   surf_lsm_h(l)%albedo(m,ind_pav_green) =                           &
                              albedo_pars(0,surf_lsm_h(l)%albedo_type(m,ind_pav_green))
                IF ( surf_lsm_h(l)%albedo_type(m,ind_wat_win) /= 0 )                 &
                   surf_lsm_h(l)%albedo(m,ind_wat_win) =                             &
                              albedo_pars(0,surf_lsm_h(l)%albedo_type(m,ind_wat_win))
             ENDDO
             DO  m = 1, surf_usm_h(l)%ns
                IF ( surf_usm_h(l)%albedo_type(m,ind_veg_wall) /= 0 )                &
                   surf_usm_h(l)%albedo(m,ind_veg_wall) =                            &
                              albedo_pars(0,surf_usm_h(l)%albedo_type(m,ind_veg_wall))
                IF ( surf_usm_h(l)%albedo_type(m,ind_pav_green) /= 0 )               &
                   surf_usm_h(l)%albedo(m,ind_pav_green) =                           &
                              albedo_pars(0,surf_usm_h(l)%albedo_type(m,ind_pav_green))
                IF ( surf_usm_h(l)%albedo_type(m,ind_wat_win) /= 0 )                 &
                   surf_usm_h(l)%albedo(m,ind_wat_win) =                             &
                              albedo_pars(0,surf_usm_h(l)%albedo_type(m,ind_wat_win))
             ENDDO
          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))
             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  l = 0, 1
                DO  m = 1, surf_lsm_h(l)%ns
                   i = surf_lsm_h(l)%i(m)
                   j = surf_lsm_h(l)%j(m)
                   IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
                      surf_lsm_h(l)%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
                      surf_lsm_h(l)%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
                      surf_lsm_h(l)%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
                   ENDIF
                ENDDO
                DO  m = 1, surf_usm_h(l)%ns
                   i = surf_usm_h(l)%i(m)
                   j = surf_usm_h(l)%j(m)
                   IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
                      surf_usm_h(l)%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
                      surf_usm_h(l)%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
                      surf_usm_h(l)%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
                   ENDIF
                ENDDO
             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 l = 0, 1
                DO  m = 1, surf_usm_h(l)%ns
                   i = surf_usm_h(l)%i(m)
                   j = surf_usm_h(l)%j(m)
                   k = surf_usm_h(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(4,is) == -surf_usm_h(l)%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(l)%albedo(m,ind_veg_wall) =                           &
                                     building_surface_pars_f%pars(ind_s_alb_b_wall,is)
                            surf_usm_h(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_h(l)%albedo(m,ind_wat_win) =                            &
                                     building_surface_pars_f%pars(ind_s_alb_b_win,is)
                            surf_usm_h(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_h(l)%albedo(m,ind_pav_green) =                          &
                                     building_surface_pars_f%pars(ind_s_alb_b_green,is)
                            surf_usm_h(l)%albedo_type(m,ind_pav_green) = 0
                         ENDIF

                         EXIT ! surface was found and processed
                      ENDIF
                   ENDDO
                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).
          DO  l = 0, 1
             ALLOCATE ( surf_lsm_h(l)%aldif(1:surf_lsm_h(l)%ns,0:2)       )
             ALLOCATE ( surf_lsm_h(l)%aldir(1:surf_lsm_h(l)%ns,0:2)       )
             ALLOCATE ( surf_lsm_h(l)%asdif(1:surf_lsm_h(l)%ns,0:2)       )
             ALLOCATE ( surf_lsm_h(l)%asdir(1:surf_lsm_h(l)%ns,0:2)       )
             ALLOCATE ( surf_lsm_h(l)%rrtm_aldif(1:surf_lsm_h(l)%ns,0:2)  )
             ALLOCATE ( surf_lsm_h(l)%rrtm_aldir(1:surf_lsm_h(l)%ns,0:2)  )
             ALLOCATE ( surf_lsm_h(l)%rrtm_asdif(1:surf_lsm_h(l)%ns,0:2)  )
             ALLOCATE ( surf_lsm_h(l)%rrtm_asdir(1:surf_lsm_h(l)%ns,0:2)  )

             ALLOCATE ( surf_usm_h(l)%aldif(1:surf_usm_h(l)%ns,0:2)       )
             ALLOCATE ( surf_usm_h(l)%aldir(1:surf_usm_h(l)%ns,0:2)       )
             ALLOCATE ( surf_usm_h(l)%asdif(1:surf_usm_h(l)%ns,0:2)       )
             ALLOCATE ( surf_usm_h(l)%asdir(1:surf_usm_h(l)%ns,0:2)       )
             ALLOCATE ( surf_usm_h(l)%rrtm_aldif(1:surf_usm_h(l)%ns,0:2)  )
             ALLOCATE ( surf_usm_h(l)%rrtm_aldir(1:surf_usm_h(l)%ns,0:2)  )
             ALLOCATE ( surf_usm_h(l)%rrtm_asdif(1:surf_usm_h(l)%ns,0:2)  )
             ALLOCATE ( surf_usm_h(l)%rrtm_asdir(1:surf_usm_h(l)%ns,0:2)  )

!
!--          Allocate broadband albedo (temporary for the current radiation
!--          implementations)
             IF ( .NOT. ALLOCATED(surf_lsm_h(l)%albedo) )                            &
                ALLOCATE( surf_lsm_h(l)%albedo(1:surf_lsm_h(l)%ns,0:2)     )
             IF ( .NOT. ALLOCATED(surf_usm_h(l)%albedo) )                            &
                ALLOCATE( surf_usm_h(l)%albedo(1:surf_usm_h(l)%ns,0:2)     )
          ENDDO
!
!--       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.
          DO  l = 0, 1
             IF ( surf_lsm_h(l)%ns > 0 )  THEN
                surf_lsm_h(l)%aldif  = albedo_lw_dif
                surf_lsm_h(l)%aldir  = albedo_lw_dir
                surf_lsm_h(l)%asdif  = albedo_sw_dif
                surf_lsm_h(l)%asdir  = albedo_sw_dir
                surf_lsm_h(l)%albedo = albedo_sw_dif
             ENDIF
             IF ( surf_usm_h(l)%ns > 0 )  THEN
                IF ( surf_usm_h(l)%albedo_from_ascii )  THEN
                   surf_usm_h(l)%aldif  = surf_usm_h(l)%albedo
                   surf_usm_h(l)%aldir  = surf_usm_h(l)%albedo
                   surf_usm_h(l)%asdif  = surf_usm_h(l)%albedo
                   surf_usm_h(l)%asdir  = surf_usm_h(l)%albedo
                ELSE
                   surf_usm_h(l)%aldif  = albedo_lw_dif
                   surf_usm_h(l)%aldir  = albedo_lw_dir
                   surf_usm_h(l)%asdif  = albedo_sw_dif
                   surf_usm_h(l)%asdir  = albedo_sw_dir
                   surf_usm_h(l)%albedo = albedo_sw_dif
                ENDIF
             ENDIF
          ENDDO
          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

!
!--       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 l = 0, 1
             DO  m = 1, surf_lsm_h(l)%ns
!
!--             Spectral albedos for vegetation/pavement/water surfaces
                DO  ind_type = 0, 2
                   IF ( surf_lsm_h(l)%albedo_type(m,ind_type) /= 0 )  THEN
                      surf_lsm_h(l)%aldif(m,ind_type) =                              &
                                  albedo_pars(1,surf_lsm_h(l)%albedo_type(m,ind_type))
                      surf_lsm_h(l)%asdif(m,ind_type) =                              &
                                  albedo_pars(2,surf_lsm_h(l)%albedo_type(m,ind_type))
                      surf_lsm_h(l)%aldir(m,ind_type) =                              &
                                  albedo_pars(1,surf_lsm_h(l)%albedo_type(m,ind_type))
                      surf_lsm_h(l)%asdir(m,ind_type) =                              &
                                  albedo_pars(2,surf_lsm_h(l)%albedo_type(m,ind_type))
                      surf_lsm_h(l)%albedo(m,ind_type) =                             &
                                  albedo_pars(0,surf_lsm_h(l)%albedo_type(m,ind_type))
                   ENDIF
                ENDDO

             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(l)%albedo_from_ascii )  THEN
                DO  m = 1, surf_usm_h(l)%ns
!
!--                Spectral albedos for wall/green/window surfaces
                   DO  ind_type = 0, 2
                      IF ( surf_usm_h(l)%albedo_type(m,ind_type) /= 0 )  THEN
                         surf_usm_h(l)%aldif(m,ind_type) =                           &
                                  albedo_pars(1,surf_usm_h(l)%albedo_type(m,ind_type))
                         surf_usm_h(l)%asdif(m,ind_type) =                           &
                                  albedo_pars(2,surf_usm_h(l)%albedo_type(m,ind_type))
                         surf_usm_h(l)%aldir(m,ind_type) =                           &
                                  albedo_pars(1,surf_usm_h(l)%albedo_type(m,ind_type))
                         surf_usm_h(l)%asdir(m,ind_type) =                           &
                                  albedo_pars(2,surf_usm_h(l)%albedo_type(m,ind_type))
                         surf_usm_h(l)%albedo(m,ind_type) =                          &
                                  albedo_pars(0,surf_usm_h(l)%albedo_type(m,ind_type))
                      ENDIF
                   ENDDO
                ENDDO
             ENDIF
          ENDDO

          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 been already 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
!
!--       Level 3 initialization at grid points where albedo type is zero.
!--       This case, spectral albedos are taken from file if available
          IF ( albedo_pars_f%from_file )  THEN
!
!--          Horizontal
             DO l = 0, 1
                DO  m = 1, surf_lsm_h(l)%ns
                   i = surf_lsm_h(l)%i(m)
                   j = surf_lsm_h(l)%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(l)%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(l)%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(l)%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(l)%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(l)%asdif(m,ind_type) =                           &
                                                albedo_pars_f%pars_xy(2,j,i)
                   ENDDO
                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(l)%albedo_from_ascii )  THEN
                   DO  m = 1, surf_usm_h(l)%ns
                      i = surf_usm_h(l)%i(m)
                      j = surf_usm_h(l)%j(m)
!
!--                   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(l)%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(l)%aldir(m,ind_veg_wall) =                       &
                                                   albedo_pars_f%pars_xy(1,j,i)
                         surf_usm_h(l)%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(l)%asdir(m,ind_veg_wall) =                       &
                                                   albedo_pars_f%pars_xy(2,j,i)
                         surf_usm_h(l)%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(l)%aldir(m,ind_pav_green) =                      &
                                                   albedo_pars_f%pars_xy(3,j,i)
                         surf_usm_h(l)%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(l)%asdir(m,ind_pav_green) =                      &
                                                   albedo_pars_f%pars_xy(4,j,i)
                         surf_usm_h(l)%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(l)%aldir(m,ind_wat_win) =                        &
                                                   albedo_pars_f%pars_xy(5,j,i)
                         surf_usm_h(l)%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(l)%asdir(m,ind_wat_win) =                        &
                                                   albedo_pars_f%pars_xy(6,j,i)
                         surf_usm_h(l)%asdif(m,ind_wat_win) =                        &
                                                   albedo_pars_f%pars_xy(6,j,i)
                      ENDIF

                   ENDDO
                ENDIF
             ENDDO
!
!--          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 been already 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 a albedo_type
!--       to zero in order to take effect.
          IF ( building_surface_pars_f%from_file )  THEN
             DO l = 0, 1
                DO  m = 1, surf_usm_h(l)%ns
                   i = surf_usm_h(l)%i(m)
                   j = surf_usm_h(l)%j(m)
                   k = surf_usm_h(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(4,is) == -surf_usm_h(l)%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(l)%albedo(m,ind_veg_wall) =                         &
                                     building_surface_pars_f%pars(ind_s_alb_b_wall,is)
                            surf_usm_h(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_h(l)%aldir(m,ind_veg_wall) =                          &
                                     building_surface_pars_f%pars(ind_s_alb_l_wall,is)
                            surf_usm_h(l)%aldif(m,ind_veg_wall) =                          &
                                     building_surface_pars_f%pars(ind_s_alb_l_wall,is)
                            surf_usm_h(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_h(l)%asdir(m,ind_veg_wall) =                          &
                                     building_surface_pars_f%pars(ind_s_alb_s_wall,is)
                            surf_usm_h(l)%asdif(m,ind_veg_wall) =                          &
                                     building_surface_pars_f%pars(ind_s_alb_s_wall,is)
                            surf_usm_h(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_h(l)%albedo(m,ind_wat_win) =                          &
                                     building_surface_pars_f%pars(ind_s_alb_b_win,is)
                            surf_usm_h(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_h(l)%aldir(m,ind_wat_win) =                           &
                                     building_surface_pars_f%pars(ind_s_alb_l_win,is)
                            surf_usm_h(l)%aldif(m,ind_wat_win) =                           &
                                     building_surface_pars_f%pars(ind_s_alb_l_win,is)
                            surf_usm_h(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_h(l)%asdir(m,ind_wat_win) =                           &
                                     building_surface_pars_f%pars(ind_s_alb_s_win,is)
                            surf_usm_h(l)%asdif(m,ind_wat_win) =                           &
                                     building_surface_pars_f%pars(ind_s_alb_s_win,is)
                            surf_usm_h(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_h(l)%albedo(m,ind_pav_green) =                        &
                                     building_surface_pars_f%pars(ind_s_alb_b_green,is)
                            surf_usm_h(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_h(l)%aldir(m,ind_pav_green) =                         &
                                     building_surface_pars_f%pars(ind_s_alb_l_green,is)
                            surf_usm_h(l)%aldif(m,ind_pav_green) =                         &
                                     building_surface_pars_f%pars(ind_s_alb_l_green,is)
                            surf_usm_h(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_h(l)%asdir(m,ind_pav_green) =                         &
                                     building_surface_pars_f%pars(ind_s_alb_s_green,is)
                            surf_usm_h(l)%asdif(m,ind_pav_green) =                         &
                                     building_surface_pars_f%pars(ind_s_alb_s_green,is)
                            surf_usm_h(l)%albedo_type(m,ind_pav_green) = 0
                         ENDIF

                         EXIT ! surface was found and processed
                      ENDIF
                   ENDDO
                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_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
             ENDDO
          ENDIF

!
!--       Calculate initial values of current (cosine of) the 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
             DO l = 0, 1
                CALL calc_albedo( surf_lsm_h(l) )
                CALL calc_albedo( surf_usm_h(l) )
             ENDDO
!
!--          Vertically aligned natural and urban surfaces
             DO  l = 0, 3
                CALL calc_albedo( surf_lsm_v(l) )
                CALL calc_albedo( surf_usm_v(l) )
             ENDDO
#endif
          ELSE
!
!--          Initialize sun-inclination independent spectral albedos
!--          Horizontal surfaces
             DO  l = 0, 1
                IF ( surf_lsm_h(l)%ns > 0 )  THEN
                   surf_lsm_h(l)%rrtm_aldir = surf_lsm_h(l)%aldir
                   surf_lsm_h(l)%rrtm_asdir = surf_lsm_h(l)%asdir
                   surf_lsm_h(l)%rrtm_aldif = surf_lsm_h(l)%aldif
                   surf_lsm_h(l)%rrtm_asdif = surf_lsm_h(l)%asdif
                ENDIF
                IF ( surf_usm_h(l)%ns > 0 )  THEN
                   surf_usm_h(l)%rrtm_aldir = surf_usm_h(l)%aldir
                   surf_usm_h(l)%rrtm_asdir = surf_usm_h(l)%asdir
                   surf_usm_h(l)%rrtm_aldif = surf_usm_h(l)%aldif
                   surf_usm_h(l)%rrtm_asdif = surf_usm_h(l)%asdif
                ENDIF
             ENDDO
!
!--          Vertical surfaces
             DO  l = 0, 3
                IF ( surf_lsm_v(l)%ns > 0 )  THEN
                   surf_lsm_v(l)%rrtm_aldir = surf_lsm_v(l)%aldir
                   surf_lsm_v(l)%rrtm_asdir = surf_lsm_v(l)%asdir
                   surf_lsm_v(l)%rrtm_aldif = surf_lsm_v(l)%aldif
                   surf_lsm_v(l)%rrtm_asdif = surf_lsm_v(l)%asdif
                ENDIF
                IF ( surf_usm_v(l)%ns > 0 )  THEN
                   surf_usm_v(l)%rrtm_aldir = surf_usm_v(l)%aldir
                   surf_usm_v(l)%rrtm_asdir = surf_usm_v(l)%asdir
                   surf_usm_v(l)%rrtm_aldif = surf_usm_v(l)%aldif
                   surf_usm_v(l)%rrtm_asdif = surf_usm_v(l)%asdif
                ENDIF
             ENDDO

          ENDIF

!
!--       Allocate 3d arrays of radiative fluxes and heating rates
          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
             ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
             rad_sw_in = 0.0_wp
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
             ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
             ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
             rad_sw_out = 0.0_wp
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
             ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_sw_hr ) )  THEN
             ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
             rad_sw_hr = 0.0_wp
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) )  THEN
             ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
             rad_sw_hr_av = 0.0_wp
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) )  THEN
             ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
             rad_sw_cs_hr = 0.0_wp
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) )  THEN
             ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
             rad_sw_cs_hr_av = 0.0_wp
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
             ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
             rad_lw_in = 0.0_wp
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
             ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
             ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
            rad_lw_out = 0.0_wp
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
             ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_lw_hr ) )  THEN
             ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
             rad_lw_hr = 0.0_wp
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) )  THEN
             ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
             rad_lw_hr_av = 0.0_wp
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) )  THEN
             ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
             rad_lw_cs_hr = 0.0_wp
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) )  THEN
             ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
             rad_lw_cs_hr_av = 0.0_wp
          ENDIF

          ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
          ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
          rad_sw_cs_in  = 0.0_wp
          rad_sw_cs_out = 0.0_wp

          ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
          ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
          rad_lw_cs_in  = 0.0_wp
          rad_lw_cs_out = 0.0_wp

!
!--       Allocate 1-element array for surface temperature
!--       (RRTMG anticipates an array as passed argument).
          ALLOCATE ( rrtm_tsfc(1) )
!
!--       Allocate surface emissivity.
!--       Values will be given directly before calling rrtm_lw.
          ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )

!
!--       Initialize RRTMG, before check if files are existent
          INQUIRE( FILE='rrtmg_lw.nc', EXIST=lw_exists )
          IF ( .NOT. lw_exists )  THEN
             message_string = 'Input file rrtmg_lw.nc' //                &
                            '&for rrtmg missing. ' // &
                            '&Please provide <jobname>_lsw file in the INPUT directory.'
             CALL message( 'radiation_init', 'PA0583', 1, 2, 0, 6, 0 )
          ENDIF
          INQUIRE( FILE='rrtmg_sw.nc', EXIST=sw_exists )
          IF ( .NOT. sw_exists )  THEN
             message_string = 'Input file rrtmg_sw.nc' //                &
                            '&for rrtmg missing. ' // &
                            '&Please provide <jobname>_rsw file in the INPUT directory.'
             CALL message( 'radiation_init', 'PA0584', 1, 2, 0, 6, 0 )
          ENDIF

          IF ( lw_radiation )  CALL rrtmg_lw_ini ( c_p )
          IF ( sw_radiation )  CALL rrtmg_sw_ini ( c_p )

!
!--       Set input files for RRTMG
          INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists)
          IF ( .NOT. snd_exists )  THEN
             rrtm_input_file = "rrtmg_lw.nc"
          ENDIF

!
!--       Read vertical layers for RRTMG from sounding data
!--       The routine provides nzt_rad, hyp_snd(1:nzt_rad),
!--       t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
!--       rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
          CALL read_sounding_data

!
!--       Read trace gas profiles from file. This routine provides
!--       the rrtm_ arrays (1:nzt_rad+1)
          CALL read_trace_gas_data
#endif
       ENDIF
!
!--    Initializaion actions exclusively required for external
!--    radiation forcing
       IF ( radiation_scheme == 'external' )  THEN
!
!--       Open the radiation input file. Note, for child domain, a dynamic
!--       input file is often not provided. In order to do not need to
!--       duplicate the dynamic input file just for the radiation input, take
!--       it from the dynamic file for the parent if not available for the
!--       child domain(s). In this case this is possible because radiation
!--       input should be the same for each model.
          INQUIRE( FILE = TRIM( input_file_dynamic ),                          &
                   EXIST = radiation_input_root_domain  )

          IF ( .NOT. input_pids_dynamic  .AND.                                 &
               .NOT. radiation_input_root_domain )  THEN
             message_string = 'In case of external radiation forcing ' //      &
                              'a dynamic input file is required. If no ' //    &
                              'dynamic input for the child domain(s) is ' //   &
                              'provided, at least one for the root domain ' // &
                              'is needed.'
             CALL message( 'radiation_init', 'PA0315', 1, 2, 0, 6, 0 )
          ENDIF
#if defined( __netcdf )
!
!--       Open dynamic input file for child domain if available, else, open
!--       dynamic input file for the root domain.
          IF ( input_pids_dynamic )  THEN
             CALL open_read_file( TRIM( input_file_dynamic ) //                &
                                  TRIM( coupling_char ),                       &
                                  pids_id )
          ELSEIF ( radiation_input_root_domain )  THEN
             CALL open_read_file( TRIM( input_file_dynamic ),                  &
                                  pids_id )
          ENDIF

          CALL inquire_num_variables( pids_id, num_var_pids )
!
!--       Allocate memory to store variable names and read them
          ALLOCATE( vars_pids(1:num_var_pids) )
          CALL inquire_variable_names( pids_id, vars_pids )
!
!--       Input time dimension.
          IF ( check_existence( vars_pids, 'time_rad' ) )  THEN
             CALL get_dimension_length( pids_id, ntime, 'time_rad' )

             ALLOCATE( time_rad_f%var1d(0:ntime-1) )
!
!--          Read variable
             CALL get_variable( pids_id, 'time_rad', time_rad_f%var1d )

             time_rad_f%from_file = .TRUE.
          ENDIF
!
!--       Input shortwave downwelling.
          IF ( check_existence( vars_pids, 'rad_sw_in' ) )  THEN
!
!--          Get _FillValue attribute
             CALL get_attribute( pids_id, char_fill, rad_sw_in_f%fill,         &
                                 .FALSE., 'rad_sw_in' )
!
!--          Get level-of-detail
             CALL get_attribute( pids_id, char_lod, rad_sw_in_f%lod,           &
                                 .FALSE., 'rad_sw_in' )
!
!--          Level-of-detail 1 - radiation depends only on time_rad
             IF ( rad_sw_in_f%lod == 1 )  THEN
                ALLOCATE( rad_sw_in_f%var1d(0:ntime-1) )
                CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var1d )
                rad_sw_in_f%from_file = .TRUE.
!
!--          Level-of-detail 2 - radiation depends on time_rad, y, x
             ELSEIF ( rad_sw_in_f%lod == 2 )  THEN
                ALLOCATE( rad_sw_in_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )

                CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var3d,    &
                                   nxl, nxr, nys, nyn, 0, ntime-1 )

                rad_sw_in_f%from_file = .TRUE.
             ELSE
                message_string = '"rad_sw_in" has no valid lod attribute'
                CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
             ENDIF
          ENDIF
!
!--       Input longwave downwelling.
          IF ( check_existence( vars_pids, 'rad_lw_in' ) )  THEN
!
!--          Get _FillValue attribute
             CALL get_attribute( pids_id, char_fill, rad_lw_in_f%fill,         &
                                 .FALSE., 'rad_lw_in' )
!
!--          Get level-of-detail
             CALL get_attribute( pids_id, char_lod, rad_lw_in_f%lod,           &
                                 .FALSE., 'rad_lw_in' )
!
!--          Level-of-detail 1 - radiation depends only on time_rad
             IF ( rad_lw_in_f%lod == 1 )  THEN
                ALLOCATE( rad_lw_in_f%var1d(0:ntime-1) )
                CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var1d )
                rad_lw_in_f%from_file = .TRUE.
!
!--          Level-of-detail 2 - radiation depends on time_rad, y, x
             ELSEIF ( rad_lw_in_f%lod == 2 )  THEN
                ALLOCATE( rad_lw_in_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )

                CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var3d,    &
                                   nxl, nxr, nys, nyn, 0, ntime-1 )

                rad_lw_in_f%from_file = .TRUE.
             ELSE
                message_string = '"rad_lw_in" has no valid lod attribute'
                CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
             ENDIF
          ENDIF
!
!--       Input shortwave downwelling, diffuse part.
          IF ( check_existence( vars_pids, 'rad_sw_in_dif' ) )  THEN
!
!--          Read _FillValue attribute
             CALL get_attribute( pids_id, char_fill, rad_sw_in_dif_f%fill,     &
                                 .FALSE., 'rad_sw_in_dif' )
!
!--          Get level-of-detail
             CALL get_attribute( pids_id, char_lod, rad_sw_in_dif_f%lod,       &
                                 .FALSE., 'rad_sw_in_dif' )
!
!--          Level-of-detail 1 - radiation depends only on time_rad
             IF ( rad_sw_in_dif_f%lod == 1 )  THEN
                ALLOCATE( rad_sw_in_dif_f%var1d(0:ntime-1) )
                CALL get_variable( pids_id, 'rad_sw_in_dif',                   &
                                   rad_sw_in_dif_f%var1d )
                rad_sw_in_dif_f%from_file = .TRUE.
!
!--          Level-of-detail 2 - radiation depends on time_rad, y, x
             ELSEIF ( rad_sw_in_dif_f%lod == 2 )  THEN
                ALLOCATE( rad_sw_in_dif_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )

                CALL get_variable( pids_id, 'rad_sw_in_dif',                   &
                                   rad_sw_in_dif_f%var3d,                      &
                                   nxl, nxr, nys, nyn, 0, ntime-1 )

                rad_sw_in_dif_f%from_file = .TRUE.
             ELSE
                message_string = '"rad_sw_in_dif" has no valid lod attribute'
                CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
             ENDIF
          ENDIF
!
!--       Finally, close the input file and deallocate temporary arrays
          DEALLOCATE( vars_pids )

          CALL close_input_file( pids_id )
#endif
!
!--       Make some consistency checks.
          IF ( .NOT. rad_sw_in_f%from_file  .OR.                               &
               .NOT. rad_lw_in_f%from_file )  THEN
             message_string = 'In case of external radiation forcing ' //      &
                              'both, rad_sw_in and rad_lw_in are required.'
             CALL message( 'radiation_init', 'PA0195', 1, 2, 0, 6, 0 )
          ENDIF

          IF ( .NOT. time_rad_f%from_file )  THEN
             message_string = 'In case of external radiation forcing ' //      &
                              'dimension time_rad is required.'
             CALL message( 'radiation_init', 'PA0196', 1, 2, 0, 6, 0 )
          ENDIF

          CALL get_date_time( 0.0_wp, second_of_day=second_of_day )

          IF ( end_time - spinup_time > time_rad_f%var1d(ntime-1) )  THEN
             message_string = 'External radiation forcing does not cover ' //  &
                              'the entire simulation time.'
             CALL message( 'radiation_init', 'PA0314', 1, 2, 0, 6, 0 )
          ENDIF
!
!--       Check for fill values in radiation
          IF ( ALLOCATED( rad_sw_in_f%var1d ) )  THEN
             IF ( ANY( rad_sw_in_f%var1d == rad_sw_in_f%fill ) )  THEN
                message_string = 'External radiation array "rad_sw_in" ' //    &
                                 'must not contain any fill values.'
                CALL message( 'radiation_init', 'PA0197', 1, 2, 0, 6, 0 )
             ENDIF
          ENDIF

          IF ( ALLOCATED( rad_lw_in_f%var1d ) )  THEN
             IF ( ANY( rad_lw_in_f%var1d == rad_lw_in_f%fill ) )  THEN
                message_string = 'External radiation array "rad_lw_in" ' //    &
                                 'must not contain any fill values.'
                CALL message( 'radiation_init', 'PA0198', 1, 2, 0, 6, 0 )
             ENDIF
          ENDIF

          IF ( ALLOCATED( rad_sw_in_dif_f%var1d ) )  THEN
             IF ( ANY( rad_sw_in_dif_f%var1d == rad_sw_in_dif_f%fill ) )  THEN
                message_string = 'External radiation array "rad_sw_in_dif" ' //&
                                 'must not contain any fill values.'
                CALL message( 'radiation_init', 'PA0199', 1, 2, 0, 6, 0 )
             ENDIF
          ENDIF

          IF ( ALLOCATED( rad_sw_in_f%var3d ) )  THEN
             IF ( ANY( rad_sw_in_f%var3d == rad_sw_in_f%fill ) )  THEN
                message_string = 'External radiation array "rad_sw_in" ' //    &
                                 'must not contain any fill values.'
                CALL message( 'radiation_init', 'PA0197', 1, 2, 0, 6, 0 )
             ENDIF
          ENDIF

          IF ( ALLOCATED( rad_lw_in_f%var3d ) )  THEN
             IF ( ANY( rad_lw_in_f%var3d == rad_lw_in_f%fill ) )  THEN
                message_string = 'External radiation array "rad_lw_in" ' //    &
                                 'must not contain any fill values.'
                CALL message( 'radiation_init', 'PA0198', 1, 2, 0, 6, 0 )
             ENDIF
          ENDIF

          IF ( ALLOCATED( rad_sw_in_dif_f%var3d ) )  THEN
             IF ( ANY( rad_sw_in_dif_f%var3d == rad_sw_in_dif_f%fill ) )  THEN
                message_string = 'External radiation array "rad_sw_in_dif" ' //&
                                 'must not contain any fill values.'
                CALL message( 'radiation_init', 'PA0199', 1, 2, 0, 6, 0 )
             ENDIF
          ENDIF
!
!--       Currently, 2D external radiation input is not possible in
!--       combination with topography where average radiation is used.
          IF ( ( rad_lw_in_f%lod == 2  .OR.  rad_sw_in_f%lod == 2  .OR.      &
                 rad_sw_in_dif_f%lod == 2  )  .AND. average_radiation )  THEN
             message_string = 'External radiation with lod = 2 is currently '//&
                              'not possible with average_radiation = .T..'
                CALL message( 'radiation_init', 'PA0670', 1, 2, 0, 6, 0 )
          ENDIF
!
!--       All radiation input should have the same level of detail. The sum
!--       of lods divided by the number of available radiation arrays must be
!--       1 (if all are lod = 1) or 2 (if all are lod = 2).
          IF ( REAL( MERGE( rad_lw_in_f%lod, 0, rad_lw_in_f%from_file ) +       &
                     MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) +       &
                     MERGE( rad_sw_in_dif_f%lod, 0, rad_sw_in_dif_f%from_file ),&
                     KIND = wp ) /                                              &
                   ( MERGE( 1.0_wp, 0.0_wp, rad_lw_in_f%from_file ) +           &
                     MERGE( 1.0_wp, 0.0_wp, rad_sw_in_f%from_file ) +           &
                     MERGE( 1.0_wp, 0.0_wp, rad_sw_in_dif_f%from_file ) )       &
                     /= 1.0_wp  .AND.                                           &
               REAL( MERGE( rad_lw_in_f%lod, 0, rad_lw_in_f%from_file ) +       &
                     MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) +       &
                     MERGE( rad_sw_in_dif_f%lod, 0, rad_sw_in_dif_f%from_file ),&
                     KIND = wp ) /                                              &
                   ( MERGE( 1.0_wp, 0.0_wp, rad_lw_in_f%from_file ) +           &
                     MERGE( 1.0_wp, 0.0_wp, rad_sw_in_f%from_file ) +           &
                     MERGE( 1.0_wp, 0.0_wp, rad_sw_in_dif_f%from_file ) )       &
                     /= 2.0_wp )  THEN
             message_string = 'External radiation input should have the same '//&
                              'lod.'
             CALL message( 'radiation_init', 'PA0673', 1, 2, 0, 6, 0 )
          ENDIF

       ENDIF
!
!--    Perform user actions if required
       CALL user_init_radiation

!
!--    Calculate radiative fluxes at model start
       SELECT CASE ( TRIM( radiation_scheme ) )

          CASE ( 'rrtmg' )
             CALL radiation_rrtmg

          CASE ( 'clear-sky' )
             CALL radiation_clearsky

          CASE ( 'constant' )
             CALL radiation_constant

          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

!
!--    Find all discretized apparent solar positions for radiation interaction.
       IF ( radiation_interactions )  CALL radiation_presimulate_solar_pos

!
!--    If required, read or calculate and write out the SVF
       IF ( radiation_interactions .AND. read_svf)  THEN
!
!--       Read sky-view factors and further required data from file
          CALL radiation_read_svf()

       ELSEIF ( radiation_interactions .AND. .NOT. read_svf)  THEN
!
!--       calculate SFV and CSF
          CALL radiation_calc_svf()
       ENDIF

       IF ( radiation_interactions .AND. write_svf)  THEN
!
!--       Write svf, csf svfsurf and csfsurf data to file
          CALL radiation_write_svf()
       ENDIF

!
!--    Adjust radiative fluxes. In case of urban and land surfaces, also
!--    call an initial interaction.
       IF ( radiation_interactions )  THEN
          CALL radiation_interaction
       ENDIF

       IF ( debug_output )  CALL debug_message( 'radiation_init', 'end' )

       RETURN !todo: remove, I don't see what we need this for here

    END SUBROUTINE radiation_init


!------------------------------------------------------------------------------!
! Description:
! ------------
!> A simple clear sky radiation model
!------------------------------------------------------------------------------!
    SUBROUTINE radiation_external

       IMPLICIT NONE

       INTEGER(iwp) ::  l   !< running index for surface orientation
       INTEGER(iwp) ::  t   !< index of current timestep
       INTEGER(iwp) ::  tm  !< index of previous timestep

       LOGICAL      ::  horizontal !< flag indicating treatment of horinzontal surfaces

       REAL(wp) ::  fac_dt               !< interpolation factor
       REAL(wp) ::  second_of_day_init   !< second of the day at model start

       TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine

!
!--    Calculate current zenith angle
       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 )
!
!--    Interpolate external radiation on current timestep
       IF ( time_since_reference_point  <= 0.0_wp )  THEN
          t      = 0
          tm     = 0
          fac_dt = 0
       ELSE
          CALL get_date_time( 0.0_wp, second_of_day=second_of_day_init )
          t = 0
          DO WHILE ( time_rad_f%var1d(t) <= time_since_reference_point )
             t = t + 1
          ENDDO

          tm = MAX( t-1, 0 )

          fac_dt = ( time_since_reference_point                                &
                   - time_rad_f%var1d(tm) + dt_3d )                            &
                 / ( time_rad_f%var1d(t)  - time_rad_f%var1d(tm) )
          fac_dt = MIN( 1.0_wp, fac_dt )
       ENDIF
!
!--    Call clear-sky calculation for each surface orientation.
!--    First, horizontal surfaces
       horizontal = .TRUE.
       DO  l = 0, 1
          surf => surf_lsm_h(l)
          CALL radiation_external_surf
          surf => surf_usm_h(l)
          CALL radiation_external_surf
       ENDDO
       horizontal = .FALSE.
!
!--    Vertical surfaces
       DO  l = 0, 3
          surf => surf_lsm_v(l)
          CALL radiation_external_surf
          surf => surf_usm_v(l)
          CALL radiation_external_surf
       ENDDO

       CONTAINS

          SUBROUTINE radiation_external_surf

             USE control_parameters

             IMPLICIT NONE

             INTEGER(iwp) ::  i    !< grid index along x-dimension
             INTEGER(iwp) ::  j    !< grid index along y-dimension
             INTEGER(iwp) ::  k    !< grid index along z-dimension
             INTEGER(iwp) ::  m    !< running index for surface elements

             REAL(wp) ::  lw_in     !< downwelling longwave radiation, interpolated value
             REAL(wp) ::  sw_in     !< downwelling shortwave radiation, interpolated value
             REAL(wp) ::  sw_in_dif !< downwelling diffuse shortwave radiation, interpolated value

             IF ( surf%ns < 1 )  RETURN
!
!--          level-of-detail = 1. Note, here it must be distinguished between
!--          averaged radiation and non-averaged radiation for the upwelling
!--          fluxes.
             IF ( rad_sw_in_f%lod == 1 )  THEN

                sw_in = ( 1.0_wp - fac_dt ) * rad_sw_in_f%var1d(tm)            &
                                   + fac_dt * rad_sw_in_f%var1d(t)

                lw_in = ( 1.0_wp - fac_dt ) * rad_lw_in_f%var1d(tm)            &
                                   + fac_dt * rad_lw_in_f%var1d(t)
!
!--             Limit shortwave incoming radiation to positive values, in order
!--             to overcome possible observation errors.
                sw_in = MAX( 0.0_wp, sw_in )
                sw_in = MERGE( sw_in, 0.0_wp, sun_up )

                surf%rad_sw_in = sw_in
                surf%rad_lw_in = lw_in

                IF ( average_radiation )  THEN
                   surf%rad_sw_out = albedo_urb * surf%rad_sw_in

                   surf%rad_lw_out = emissivity_urb * sigma_sb * t_rad_urb**4  &
                                  + ( 1.0_wp - emissivity_urb ) * surf%rad_lw_in

                   surf%rad_net = surf%rad_sw_in - surf%rad_sw_out             &
                                + surf%rad_lw_in - surf%rad_lw_out

                   surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb          &
                                                     * sigma_sb                &
                                                     * t_rad_urb**3
                ELSE
                   DO  m = 1, surf%ns
                      k = surf%k(m)
                      surf%rad_sw_out(m) = ( surf%frac(m,ind_veg_wall)  *      &
                                             surf%albedo(m,ind_veg_wall)       &
                                           + surf%frac(m,ind_pav_green) *      &
                                             surf%albedo(m,ind_pav_green)      &
                                           + surf%frac(m,ind_wat_win)   *      &
                                             surf%albedo(m,ind_wat_win) )      &
                                           * surf%rad_sw_in(m)

                      surf%rad_lw_out(m) = ( surf%frac(m,ind_veg_wall)  *      &
                                             surf%emissivity(m,ind_veg_wall)   &
                                           + surf%frac(m,ind_pav_green) *      &
                                             surf%emissivity(m,ind_pav_green)  &
                                           + surf%frac(m,ind_wat_win)   *      &
                                             surf%emissivity(m,ind_wat_win)    &
                                           )                                   &
                                           * sigma_sb                          &
                                           * ( surf%pt_surface(m) * exner(k) )**4

                      surf%rad_lw_out_change_0(m) =                            &
                                         ( surf%frac(m,ind_veg_wall)  *        &
                                           surf%emissivity(m,ind_veg_wall)     &
                                         + surf%frac(m,ind_pav_green) *        &
                                           surf%emissivity(m,ind_pav_green)    &
                                         + surf%frac(m,ind_wat_win)   *        &
                                           surf%emissivity(m,ind_wat_win)      &
                                         ) * 4.0_wp * sigma_sb                 &
                                         * ( surf%pt_surface(m) * exner(k) )**3
                   ENDDO

                ENDIF
!
!--             If diffuse shortwave radiation is available, store it on
!--             the respective files.
                IF ( rad_sw_in_dif_f%from_file )  THEN
                   sw_in_dif= ( 1.0_wp - fac_dt ) * rad_sw_in_dif_f%var1d(tm)  &
                                         + fac_dt * rad_sw_in_dif_f%var1d(t)

                   IF ( ALLOCATED( rad_sw_in_diff ) )  rad_sw_in_diff = sw_in_dif
                   IF ( ALLOCATED( rad_sw_in_dir  ) )  rad_sw_in_dir  = sw_in  &
                                                                    - sw_in_dif
!
!--                Diffuse longwave radiation equals the total downwelling
!--                longwave radiation
                   IF ( ALLOCATED( rad_lw_in_diff ) )  rad_lw_in_diff = lw_in
                ENDIF
!
!--          level-of-detail = 2
             ELSE

                DO  m = 1, surf%ns
                   i = surf%i(m)
                   j = surf%j(m)
                   k = surf%k(m)

                   surf%rad_sw_in(m) = ( 1.0_wp - fac_dt )                     &
                                            * rad_sw_in_f%var3d(tm,j,i)        &
                                   + fac_dt * rad_sw_in_f%var3d(t,j,i)
!
!--                Limit shortwave incoming radiation to positive values, in
!--                order to overcome possible observation errors.
                   surf%rad_sw_in(m) = MAX( 0.0_wp, surf%rad_sw_in(m) )
                   surf%rad_sw_in(m) = MERGE( surf%rad_sw_in(m), 0.0_wp, sun_up )

                   surf%rad_lw_in(m) = ( 1.0_wp - fac_dt )                     &
                                            * rad_lw_in_f%var3d(tm,j,i)        &
                                   + fac_dt * rad_lw_in_f%var3d(t,j,i)
!
!--                Weighted average according to surface fraction.
                   surf%rad_sw_out(m) = ( surf%frac(m,ind_veg_wall)  *         &
                                          surf%albedo(m,ind_veg_wall)          &
                                        + surf%frac(m,ind_pav_green) *         &
                                          surf%albedo(m,ind_pav_green)         &
                                        + surf%frac(m,ind_wat_win)   *         &
                                          surf%albedo(m,ind_wat_win) )         &
                                        * surf%rad_sw_in(m)

                   surf%rad_lw_out(m) = ( surf%frac(m,ind_veg_wall)  *         &
                                          surf%emissivity(m,ind_veg_wall)      &
                                        + surf%frac(m,ind_pav_green) *         &
                                          surf%emissivity(m,ind_pav_green)     &
                                        + surf%frac(m,ind_wat_win)   *         &
                                          surf%emissivity(m,ind_wat_win)       &
                                        )                                      &
                                        * sigma_sb                             &
                                        * ( surf%pt_surface(m) * exner(k) )**4

                   surf%rad_lw_out_change_0(m) =                               &
                                      ( surf%frac(m,ind_veg_wall)  *           &
                                        surf%emissivity(m,ind_veg_wall)        &
                                      + surf%frac(m,ind_pav_green) *           &
                                        surf%emissivity(m,ind_pav_green)       &
                                      + surf%frac(m,ind_wat_win)   *           &
                                        surf%emissivity(m,ind_wat_win)         &
                                      ) * 4.0_wp * sigma_sb                    &
                                      * ( surf%pt_surface(m) * exner(k) )**3

                   surf%rad_net(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m)    &
                                   + surf%rad_lw_in(m) - surf%rad_lw_out(m)
!
!--                If diffuse shortwave radiation is available, store it on
!--                the respective files.
                   IF ( rad_sw_in_dif_f%from_file )  THEN
                      IF ( ALLOCATED( rad_sw_in_diff ) )                       &
                         rad_sw_in_diff(j,i) = ( 1.0_wp - fac_dt )             &
                                              * rad_sw_in_dif_f%var3d(tm,j,i)  &
                                     + fac_dt * rad_sw_in_dif_f%var3d(t,j,i)
!
!--                   dir = sw_in - sw_in_dif.
                      IF ( ALLOCATED( rad_sw_in_dir  ) )                       &
                         rad_sw_in_dir(j,i)  = surf%rad_sw_in(m) -             &
                                               rad_sw_in_diff(j,i)
!
!--                   Diffuse longwave radiation equals the total downwelling
!--                   longwave radiation
                      IF ( ALLOCATED( rad_lw_in_diff ) )                       &
                         rad_lw_in_diff(j,i) = surf%rad_lw_in(m)
                   ENDIF

                ENDDO

             ENDIF
!
!--          Store radiation also on 2D arrays, which are still used for
!--          direct-diffuse splitting. Note, this is only required
!--          for horizontal surfaces, which covers all x,y position.
             IF ( horizontal )  THEN
                DO  m = 1, surf%ns
                   i = surf%i(m)
                   j = surf%j(m)

                   rad_sw_in(0,j,i)  = surf%rad_sw_in(m)
                   rad_lw_in(0,j,i)  = surf%rad_lw_in(m)
                   rad_sw_out(0,j,i) = surf%rad_sw_out(m)
                   rad_lw_out(0,j,i) = surf%rad_lw_out(m)
                ENDDO
             ENDIF

          END SUBROUTINE radiation_external_surf

    END SUBROUTINE radiation_external

!------------------------------------------------------------------------------!
! Description:
! ------------
!> A simple clear sky radiation model
!------------------------------------------------------------------------------!
    SUBROUTINE radiation_clearsky

       IMPLICIT NONE

       INTEGER(iwp) ::  l         !< running index for surface orientation

       LOGICAL      ::  horizontal !< flag indicating treatment of horinzontal surfaces

       REAL(wp)     ::  pt1       !< potential temperature at first grid level or mean value at urban layer top
       REAL(wp)     ::  pt1_l     !< potential temperature at first grid level or mean value at urban layer top at local subdomain
       REAL(wp)     ::  ql1       !< liquid water mixing ratio at first grid level or mean value at urban layer top
       REAL(wp)     ::  ql1_l     !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain

       TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine

!
!--    Calculate current zenith angle
       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 sky transmissivity
       sky_trans = 0.6_wp + 0.2_wp * cos_zenith

!
!--    Calculate value of the Exner function at model surface
!
!--    In case averaged radiation is used, calculate mean temperature and
!--    liquid water mixing ratio at the urban-layer top.
       IF ( average_radiation ) THEN
          pt1   = 0.0_wp
          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1   = 0.0_wp

          pt1_l = SUM( pt(nz_urban_t,nys:nyn,nxl:nxr) )
          IF ( bulk_cloud_model  .OR.  cloud_droplets  )  ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) )

#if defined( __parallel )
          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
          CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
          IF ( ierr /= 0 ) THEN
              WRITE(9,*) 'Error MPI_AllReduce1:', ierr, pt1_l, pt1
              FLUSH(9)
          ENDIF

          IF ( bulk_cloud_model  .OR.  cloud_droplets ) THEN
              CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
              IF ( ierr /= 0 ) THEN
                  WRITE(9,*) 'Error MPI_AllReduce2:', ierr, ql1_l, ql1
                  FLUSH(9)
              ENDIF
          ENDIF
#else
          pt1 = pt1_l
          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1 = ql1_l
#endif

          IF ( bulk_cloud_model  .OR.  cloud_droplets  )  pt1 = pt1 + lv_d_cp / exner(nz_urban_t) * ql1
!
!--       Finally, divide by number of grid points
          pt1 = pt1 / REAL( ( nx + 1 ) * ( ny + 1 ), KIND=wp )
       ENDIF
!
!--    Call clear-sky calculation for each surface orientation.
!--    First, horizontal surfaces
       horizontal = .TRUE.
       DO  l = 0, 1
          surf => surf_lsm_h(l)
          CALL radiation_clearsky_surf
          surf => surf_usm_h(l)
          CALL radiation_clearsky_surf
       ENDDO
       horizontal = .FALSE.
!
!--    Vertical surfaces
       DO  l = 0, 3
          surf => surf_lsm_v(l)
          CALL radiation_clearsky_surf
          surf => surf_usm_v(l)
          CALL radiation_clearsky_surf
       ENDDO

       CONTAINS

          SUBROUTINE radiation_clearsky_surf

             IMPLICIT NONE

             INTEGER(iwp) ::  i         !< index x-direction
             INTEGER(iwp) ::  j         !< index y-direction
             INTEGER(iwp) ::  k         !< index z-direction
             INTEGER(iwp) ::  m         !< running index for surface elements

             IF ( surf%ns < 1 )  RETURN

!
!--          Calculate radiation fluxes and net radiation (rad_net) assuming
!--          homogeneous urban radiation conditions.
             IF ( average_radiation ) THEN

                k = nz_urban_t

                surf%rad_sw_in  = solar_constant * sky_trans * cos_zenith
                surf%rad_sw_out = albedo_urb * surf%rad_sw_in

                surf%rad_lw_in  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k+1))**4

                surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4   &
                                    + (1.0_wp - emissivity_urb) * surf%rad_lw_in

                surf%rad_net = surf%rad_sw_in - surf%rad_sw_out                &
                             + surf%rad_lw_in - surf%rad_lw_out

                surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb * sigma_sb  &
                                           * (t_rad_urb)**3

!
!--          Calculate radiation fluxes and net radiation (rad_net) for each surface
!--          element.
             ELSE

                DO  m = 1, surf%ns
                   i = surf%i(m)
                   j = surf%j(m)
                   k = surf%k(m)

                   surf%rad_sw_in(m) = solar_constant * sky_trans * cos_zenith

!
!--                Weighted average according to surface fraction.
!--                ATTENTION: when radiation interactions are switched on the
!--                calculated fluxes below are not actually used as they are
!--                overwritten in radiation_interaction.
                   surf%rad_sw_out(m) = ( surf%frac(m,ind_veg_wall)  *         &
                                          surf%albedo(m,ind_veg_wall)          &
                                        + surf%frac(m,ind_pav_green) *         &
                                          surf%albedo(m,ind_pav_green)         &
                                        + surf%frac(m,ind_wat_win)   *         &
                                          surf%albedo(m,ind_wat_win) )         &
                                        * surf%rad_sw_in(m)

                   surf%rad_lw_out(m) = ( surf%frac(m,ind_veg_wall)  *         &
                                          surf%emissivity(m,ind_veg_wall)      &
                                        + surf%frac(m,ind_pav_green) *         &
                                          surf%emissivity(m,ind_pav_green)     &
                                        + surf%frac(m,ind_wat_win)   *         &
                                          surf%emissivity(m,ind_wat_win)       &
                                        )                                      &
                                        * sigma_sb                             &
                                        * ( surf%pt_surface(m) * exner(nzb) )**4

                   surf%rad_lw_out_change_0(m) =                               &
                                      ( surf%frac(m,ind_veg_wall)  *           &
                                        surf%emissivity(m,ind_veg_wall)        &
                                      + surf%frac(m,ind_pav_green) *           &
                                        surf%emissivity(m,ind_pav_green)       &
                                      + surf%frac(m,ind_wat_win)   *           &
                                        surf%emissivity(m,ind_wat_win)         &
                                      ) * 4.0_wp * sigma_sb                    &
                                      * ( surf%pt_surface(m) * exner(nzb) )** 3


                   IF ( bulk_cloud_model  .OR.  cloud_droplets  )  THEN
                      pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i)
                      surf%rad_lw_in(m)  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k))**4
                   ELSE
                      surf%rad_lw_in(m)  = emissivity_atm_clsky * sigma_sb * (pt(k,j,i) * exner(k))**4
                   ENDIF

                   surf%rad_net(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m)    &
                                   + surf%rad_lw_in(m) - surf%rad_lw_out(m)

                ENDDO

             ENDIF

!
!--          Fill out values in radiation arrays. Note, this is only required
!--          for horizontal surfaces, which covers all x,y position.
             IF ( horizontal )  THEN
                DO  m = 1, surf%ns
                   i = surf%i(m)
                   j = surf%j(m)
                   rad_sw_in(0,j,i)  = surf%rad_sw_in(m)
                   rad_sw_out(0,j,i) = surf%rad_sw_out(m)
                   rad_lw_in(0,j,i)  = surf%rad_lw_in(m)
                   rad_lw_out(0,j,i) = surf%rad_lw_out(m)
                ENDDO
             ENDIF

          END SUBROUTINE radiation_clearsky_surf

    END SUBROUTINE radiation_clearsky


!------------------------------------------------------------------------------!
! Description:
! ------------
!> This scheme keeps the prescribed net radiation constant during the run
!------------------------------------------------------------------------------!
    SUBROUTINE radiation_constant


       IMPLICIT NONE

       INTEGER(iwp) ::  l         !< running index for surface orientation

       LOGICAL      ::  horizontal !< flag indicating treatment of horinzontal surfaces

       REAL(wp)     ::  pt1       !< potential temperature at first grid level or mean value at urban layer top
       REAL(wp)     ::  pt1_l     !< potential temperature at first grid level or mean value at urban layer top at local subdomain
       REAL(wp)     ::  ql1       !< liquid water mixing ratio at first grid level or mean value at urban layer top
       REAL(wp)     ::  ql1_l     !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain

       TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine

!
!--    In case averaged radiation is used, calculate mean temperature and
!--    liquid water mixing ratio at the urban-layer top.
       IF ( average_radiation ) THEN
          pt1   = 0.0_wp
          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1   = 0.0_wp

          pt1_l = SUM( pt(nz_urban_t,nys:nyn,nxl:nxr) )
          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) )

#if defined( __parallel )
          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
          CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
          IF ( ierr /= 0 ) THEN
              WRITE(9,*) 'Error MPI_AllReduce3:', ierr, pt1_l, pt1
              FLUSH(9)
          ENDIF
          IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
             CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
             IF ( ierr /= 0 ) THEN
                 WRITE(9,*) 'Error MPI_AllReduce4:', ierr, ql1_l, ql1
                 FLUSH(9)
             ENDIF
          ENDIF
#else
          pt1 = pt1_l
          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1 = ql1_l
#endif
          IF ( bulk_cloud_model  .OR.  cloud_droplets )  pt1 = pt1 + lv_d_cp / exner(nz_urban_t+1) * ql1
!
!--       Finally, divide by number of grid points
          pt1 = pt1 / REAL( ( nx + 1 ) * ( ny + 1 ), KIND=wp )
       ENDIF

!
!--    First, horizontal surfaces
       horizontal = .TRUE.
       DO  l = 0, 1
          surf => surf_lsm_h(l)
          CALL radiation_constant_surf
          surf => surf_usm_h(l)
          CALL radiation_constant_surf
       ENDDO
       horizontal = .FALSE.
!
!--    Vertical surfaces
       DO  l = 0, 3
          surf => surf_lsm_v(l)
          CALL radiation_constant_surf
          surf => surf_usm_v(l)
          CALL radiation_constant_surf
       ENDDO

       CONTAINS

          SUBROUTINE radiation_constant_surf

             IMPLICIT NONE

             INTEGER(iwp) ::  i         !< index x-direction
             INTEGER(iwp) ::  ioff      !< offset between surface element and adjacent grid point along x
             INTEGER(iwp) ::  j         !< index y-direction
             INTEGER(iwp) ::  joff      !< offset between surface element and adjacent grid point along y
             INTEGER(iwp) ::  k         !< index z-direction
             INTEGER(iwp) ::  koff      !< offset between surface element and adjacent grid point along z
             INTEGER(iwp) ::  m         !< running index for surface elements

             IF ( surf%ns < 1 )  RETURN

!--          Calculate homogenoeus urban radiation fluxes
             IF ( average_radiation ) THEN

                surf%rad_net = net_radiation

                surf%rad_lw_in  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(nz_urban_t+1))**4

                surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4   &
                                    + ( 1.0_wp - emissivity_urb )             & ! shouldn't be this a bulk value -- emissivity_urb?
                                    * surf%rad_lw_in

                surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb * sigma_sb  &
                                           * t_rad_urb**3

                surf%rad_sw_in = ( surf%rad_net - surf%rad_lw_in               &
                                     + surf%rad_lw_out )                       &
                                     / ( 1.0_wp - albedo_urb )

                surf%rad_sw_out =  albedo_urb * surf%rad_sw_in

!
!--          Calculate radiation fluxes for each surface element
             ELSE
!
!--             Determine index offset between surface element and adjacent
!--             atmospheric grid point
                ioff = surf%ioff
                joff = surf%joff
                koff = surf%koff

!
!--             Prescribe net radiation and estimate the remaining radiative fluxes
                DO  m = 1, surf%ns
                   i = surf%i(m)
                   j = surf%j(m)
                   k = surf%k(m)

                   surf%rad_net(m) = net_radiation

                   IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
                      pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i)
                      surf%rad_lw_in(m)  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k))**4
                   ELSE
                      surf%rad_lw_in(m)  = emissivity_atm_clsky * sigma_sb *                 &
                                             ( pt(k,j,i) * exner(k) )**4
                   ENDIF

!
!--                Weighted average according to surface fraction.
                   surf%rad_lw_out(m) = ( surf%frac(m,ind_veg_wall)  *         &
                                          surf%emissivity(m,ind_veg_wall)      &
                                        + surf%frac(m,ind_pav_green) *         &
                                          surf%emissivity(m,ind_pav_green)     &
                                        + surf%frac(m,ind_wat_win)   *         &
                                          surf%emissivity(m,ind_wat_win)       &
                                        )                                      &
                                      * sigma_sb                               &
                                      * ( surf%pt_surface(m) * exner(nzb) )**4

                   surf%rad_sw_in(m) = ( surf%rad_net(m) - surf%rad_lw_in(m)   &
                                       + surf%rad_lw_out(m) )                  &
                                       / ( 1.0_wp -                            &
                                          ( surf%frac(m,ind_veg_wall)  *       &
                                            surf%albedo(m,ind_veg_wall)        &
                                         +  surf%frac(m,ind_pav_green) *       &
                                            surf%albedo(m,ind_pav_green)       &
                                         +  surf%frac(m,ind_wat_win)   *       &
                                            surf%albedo(m,ind_wat_win) )       &
                                         )

                   surf%rad_sw_out(m) = ( surf%frac(m,ind_veg_wall)  *         &
                                          surf%albedo(m,ind_veg_wall)          &
                                        + surf%frac(m,ind_pav_green) *         &
                                          surf%albedo(m,ind_pav_green)         &
                                        + surf%frac(m,ind_wat_win)   *         &
                                          surf%albedo(m,ind_wat_win) )         &
                                      * surf%rad_sw_in(m)

                ENDDO

             ENDIF

!
!--          Fill out values in radiation arrays. Note, this is only required
!--          for horizontal surfaces, which covers all x,y position.
             IF ( horizontal )  THEN
                DO  m = 1, surf%ns
                   i = surf%i(m)
                   j = surf%j(m)
                   rad_sw_in(0,j,i)  = surf%rad_sw_in(m)
                   rad_sw_out(0,j,i) = surf%rad_sw_out(m)
                   rad_lw_in(0,j,i)  = surf%rad_lw_in(m)
                   rad_lw_out(0,j,i) = surf%rad_lw_out(m)
                ENDDO
             ENDIF

          END SUBROUTINE radiation_constant_surf


    END SUBROUTINE radiation_constant

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Header output for radiation model
!------------------------------------------------------------------------------!
    SUBROUTINE radiation_header ( io )


       IMPLICIT NONE

       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file



!
!--    Write radiation model header
       WRITE( io, 3 )

       IF ( radiation_scheme == "constant" )  THEN
          WRITE( io, 4 ) net_radiation
       ELSEIF ( radiation_scheme == "clear-sky" )  THEN
          WRITE( io, 5 )
       ELSEIF ( radiation_scheme == "rrtmg" )  THEN
          WRITE( io, 6 )
          IF ( .NOT. lw_radiation )  WRITE( io, 10 )
          IF ( .NOT. sw_radiation )  WRITE( io, 11 )
       ELSEIF ( radiation_scheme == "external" )  THEN
          WRITE( io, 14 )
       ENDIF

       IF ( albedo_type_f%from_file  .OR.  vegetation_type_f%from_file  .OR.   &
            pavement_type_f%from_file  .OR.  water_type_f%from_file  .OR.      &
            building_type_f%from_file )  THEN
             WRITE( io, 13 )
       ELSE
          IF ( albedo_type == 0 )  THEN
             WRITE( io, 7 ) albedo
          ELSE
             WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
          ENDIF
       ENDIF
       IF ( constant_albedo )  THEN
          WRITE( io, 9 )
       ENDIF

       WRITE( io, 12 ) dt_radiation


 3 FORMAT (//' Radiation model information:'/                                  &
              ' ----------------------------'/)
 4 FORMAT ('    --> Using constant net radiation: net_radiation = ', F6.2,     &
           // 'W/m**2')
 5 FORMAT ('    --> Simple radiation scheme for clear sky is used (no clouds,',&
                   ' default)')
 6 FORMAT ('    --> RRTMG scheme is used')
 7 FORMAT (/'    User-specific surface albedo: albedo =', F6.3)
 8 FORMAT (/'    Albedo is set for land surface type: ', A)
 9 FORMAT (/'    --> Albedo is fixed during the run')
10 FORMAT (/'    --> Longwave radiation is disabled')
11 FORMAT (/'    --> Shortwave radiation is disabled.')
12 FORMAT  ('    Timestep: dt_radiation = ', F6.2, '  s')
13 FORMAT (/'    Albedo is set individually for each xy-location, according ', &
                 'to given surface type.')
14 FORMAT ('    --> External radiation forcing is used')


    END SUBROUTINE radiation_header


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Parin for &radiation_parameters for radiation model and RTM
!------------------------------------------------------------------------------!
    SUBROUTINE radiation_parin


       IMPLICIT NONE

       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file

       NAMELIST /radiation_par/   albedo, albedo_lw_dif, albedo_lw_dir,         &
                                  albedo_sw_dif, albedo_sw_dir, albedo_type,    &
                                  bufsize_alltoall,                             &
                                  constant_albedo, dt_radiation, emissivity,    &
                                  lw_radiation, max_raytracing_dist,            &
                                  min_irrf_value, mrt_geom, mrt_geom_params,    &
                                  mrt_include_sw, mrt_nlevels,                  &
                                  mrt_skip_roof, net_radiation, nrefsteps,      &
                                  plant_lw_interact, rad_angular_discretization,&
                                  radiation_interactions_on, radiation_scheme,  &
                                  raytrace_discrete_azims,                      &
                                  raytrace_discrete_elevs, raytrace_mpi_rma,    &
                                  trace_fluxes_above,                           &
                                  skip_time_do_radiation, surface_reflections,  &
                                  svfnorm_report_thresh, sw_radiation,          &
                                  unscheduled_radiation_calls


       NAMELIST /radiation_parameters/ albedo, albedo_lw_dif, albedo_lw_dir,    &
                                  albedo_sw_dif, albedo_sw_dir, albedo_type,    &
                                  bufsize_alltoall,                             &
                                  constant_albedo, dt_radiation, emissivity,    &
                                  lw_radiation, max_raytracing_dist,            &
                                  min_irrf_value, mrt_geom, mrt_geom_params,    &
                                  mrt_include_sw, mrt_nlevels,                  &
                                  mrt_skip_roof, net_radiation, nrefsteps,      &
                                  plant_lw_interact, rad_angular_discretization,&
                                  radiation_interactions_on, radiation_scheme,  &
                                  raytrace_discrete_azims,                      &
                                  raytrace_discrete_elevs, raytrace_mpi_rma,    &
                                  trace_fluxes_above,                           &
                                  skip_time_do_radiation, surface_reflections,  &
                                  svfnorm_report_thresh, sw_radiation,          &
                                  unscheduled_radiation_calls

       line = ' '

!
!--    Try to find radiation model namelist
       REWIND ( 11 )
       line = ' '
       DO WHILE ( INDEX( line, '&radiation_parameters' ) == 0 )
          READ ( 11, '(A)', END=12 )  line
       ENDDO
       BACKSPACE ( 11 )

!
!--    Read user-defined namelist
       READ ( 11, radiation_parameters, ERR = 10 )

!
!--    Set flag that indicates that the radiation model is switched on
       radiation = .TRUE.

       GOTO 14

 10    BACKSPACE( 11 )
       READ( 11 , '(A)') line
       CALL parin_fail_message( 'radiation_parameters', line )
!
!--    Try to find old namelist
 12    REWIND ( 11 )
       line = ' '
       DO WHILE ( INDEX( line, '&radiation_par' ) == 0 )
          READ ( 11, '(A)', END=14 )  line
       ENDDO
       BACKSPACE ( 11 )

!
!--    Read user-defined namelist
       READ ( 11, radiation_par, ERR = 13, END = 14 )

       message_string = 'namelist radiation_par is deprecated and will be ' // &
                     'removed in near future. Please use namelist ' //         &
                     'radiation_parameters instead'
       CALL message( 'radiation_parin', 'PA0487', 0, 1, 0, 6, 0 )

!
!--    Set flag that indicates that the radiation model is switched on
       radiation = .TRUE.

       IF ( .NOT.  radiation_interactions_on  .AND.  surface_reflections )  THEN
          message_string = 'surface_reflections is allowed only when '      // &
               'radiation_interactions_on is set to TRUE'
          CALL message( 'radiation_parin', 'PA0293',1, 2, 0, 6, 0 )
       ENDIF

       GOTO 14

 13    BACKSPACE( 11 )
       READ( 11 , '(A)') line
       CALL parin_fail_message( 'radiation_par', line )

 14    CONTINUE

    END SUBROUTINE radiation_parin


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Implementation of the RRTMG radiation_scheme
!------------------------------------------------------------------------------!
    SUBROUTINE radiation_rrtmg

#if defined ( __rrtmg )
       USE exchange_horiz_mod,                                                 &
           ONLY:  exchange_horiz

       USE indices,                                                            &
           ONLY:  nbgp

       USE palm_date_time_mod,                                                 &
           ONLY:  hours_per_day

       USE particle_attributes,                                                &
           ONLY:  grid_particles, number_of_particles, particles, prt_count

       IMPLICIT NONE


       INTEGER(iwp) ::  i, j, k, l, m, n !< loop indices
       INTEGER(iwp) ::  k_topo_l   !< topography top index
       INTEGER(iwp) ::  k_topo     !< topography top index

       REAL(wp)     ::  d_hours_day  !< 1 / hours-per-day
       REAL(wp)     ::  nc_rad, &    !< number concentration of cloud droplets
                        s_r2,   &    !< weighted sum over all droplets with r^2
                        s_r3         !< weighted sum over all droplets with r^3

       REAL(wp), DIMENSION(0:nzt+1) :: pt_av, q_av, ql_av
       REAL(wp), DIMENSION(0:0)     :: zenith   !< to provide indexed array
!
!--    Just dummy arguments
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rrtm_lw_taucld_dum,          &
                                                  rrtm_lw_tauaer_dum,          &
                                                  rrtm_sw_taucld_dum,          &
                                                  rrtm_sw_ssacld_dum,          &
                                                  rrtm_sw_asmcld_dum,          &
                                                  rrtm_sw_fsfcld_dum,          &
                                                  rrtm_sw_tauaer_dum,          &
                                                  rrtm_sw_ssaaer_dum,          &
                                                  rrtm_sw_asmaer_dum,          &
                                                  rrtm_sw_ecaer_dum

!
!--    Pre-calculate parameters
       d_hours_day = 1.0_wp / REAL( hours_per_day, KIND=wp )

!
!--    Calculate 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 )
       zenith(0) = cos_zenith
!
!--    Calculate surface albedo. In case average radiation is applied,
!--    this is not required.
#if defined( __netcdf )
       IF ( .NOT. constant_albedo )  THEN
!
!--       Horizontally aligned default, natural and urban surfaces
          DO  l = 0, 1
             CALL calc_albedo( surf_lsm_h(l) )
             CALL calc_albedo( surf_usm_h(l) )
          ENDDO
!
!--       Vertically aligned default, natural and urban surfaces
          DO  l = 0, 3
             CALL calc_albedo( surf_lsm_v(l) )
             CALL calc_albedo( surf_usm_v(l) )
          ENDDO
       ENDIF
#endif

!
!--    Prepare input data for RRTMG

!
!--    In case of large scale forcing with surface data, calculate new pressure
!--    profile. nzt_rad might be modified by these calls and all required arrays
!--    will then be re-allocated
       IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
          CALL read_sounding_data
          CALL read_trace_gas_data
       ENDIF


       IF ( average_radiation ) THEN
!
!--       Determine minimum topography top index.
          k_topo_l = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
#if defined( __parallel )
          CALL MPI_ALLREDUCE( k_topo_l, k_topo, 1, MPI_INTEGER, MPI_MIN, &
                              comm2d, ierr)
#else
          k_topo = k_topo_l
#endif

          rrtm_asdir(1)  = albedo_urb
          rrtm_asdif(1)  = albedo_urb
          rrtm_aldir(1)  = albedo_urb
          rrtm_aldif(1)  = albedo_urb

          rrtm_emis = emissivity_urb
!
!--       Calculate mean pt profile.
          CALL calc_mean_profile( pt, 4, .TRUE. )
          pt_av = hom(:, 1, 4, 0)

          IF ( humidity )  THEN
             CALL calc_mean_profile( q, 41, .TRUE. )
             q_av  = hom(:, 1, 41, 0)
          ENDIF
!
!--       Prepare profiles of temperature and H2O volume mixing ratio
          rrtm_tlev(0,k_topo+1) = t_rad_urb

          IF ( bulk_cloud_model )  THEN

             CALL calc_mean_profile( ql, 54, .TRUE. )
             ! average ql is now in hom(:, 1, 54, 0)
             ql_av = hom(:, 1, 54, 0)

             DO k = nzb+1, nzt+1
                rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp       &
                                 )**.286_wp + lv_d_cp * ql_av(k)
                rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q_av(k) - ql_av(k))
             ENDDO
          ELSE
             DO k = nzb+1, nzt+1
                rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp       &
                                 )**.286_wp
             ENDDO

             IF ( humidity )  THEN
                DO k = nzb+1, nzt+1
                   rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q_av(k)
                ENDDO
             ELSE
                rrtm_h2ovmr(0,nzb+1:nzt+1) = 0.0_wp
             ENDIF
          ENDIF

!
!--       Avoid temperature/humidity jumps at the top of the PALM domain by
!--       linear interpolation from nzt+2 to nzt+7. Jumps are induced by
!--       discrepancies between the values in the  domain and those above that
!--       are prescribed in RRTMG
          DO k = nzt+2, nzt+7
             rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1)                            &
                           + ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) )    &
                           / ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )    &
                           * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )

             rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1)                        &
                           + ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
                           / ( rrtm_play(0,nzt+8)   - rrtm_play(0,nzt+1)   )&
                           * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )

          ENDDO

!--       Linear interpolate to zw grid. Loop reaches one level further up
!--       due to the staggered grid in RRTMG
          DO k = k_topo+2, nzt+8
             rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) -        &
                                rrtm_tlay(0,k-1))                           &
                                / ( rrtm_play(0,k) - rrtm_play(0,k-1) )     &
                                * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
          ENDDO
!
!--       Calculate liquid water path and cloud fraction for each column.
!--       Note that LWP is required in g/m2 instead of kg/kg m.
          rrtm_cldfr  = 0.0_wp
          rrtm_reliq  = 0.0_wp
          rrtm_cliqwp = 0.0_wp
          rrtm_icld   = 0

          IF ( bulk_cloud_model )  THEN
             DO k = nzb+1, nzt+1
                rrtm_cliqwp(0,k) =  ql_av(k) * 1000._wp *                   &
                                    (rrtm_plev(0,k) - rrtm_plev(0,k+1))     &
                                    * 100._wp / g

                IF ( rrtm_cliqwp(0,k) > 0._wp )  THEN
                   rrtm_cldfr(0,k) = 1._wp
                   IF ( rrtm_icld == 0 )  rrtm_icld = 1

!
!--                Calculate cloud droplet effective radius
                   rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql_av(k)         &
                                     * rho_surface                          &
                                     / ( 4.0_wp * pi * nc_const * rho_l )   &
                                     )**0.33333333333333_wp                 &
                                     * EXP( LOG( sigma_gc )**2 )
!
!--                Limit effective radius
                   IF ( rrtm_reliq(0,k) > 0.0_wp )  THEN
                      rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
                      rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
                   ENDIF
                ENDIF
             ENDDO
          ENDIF

!
!--       Set surface temperature
          rrtm_tsfc = t_rad_urb

          IF ( lw_radiation )  THEN
!
!--          Due to technical reasons, copy optical depth to dummy arguments
!--          which are allocated on the exact size as the rrtmg_lw is called.
!--          As one dimesion is allocated with zero size, compiler complains
!--          that rank of the array does not match that of the
!--          assumed-shaped arguments in the RRTMG library. In order to
!--          avoid this, write to dummy arguments and give pass the entire
!--          dummy array. Seems to be the only existing work-around.
             ALLOCATE( rrtm_lw_taucld_dum(1:nbndlw+1,0:0,k_topo+1:nzt_rad+1) )
             ALLOCATE( rrtm_lw_tauaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1) )

             rrtm_lw_taucld_dum =                                              &
                             rrtm_lw_taucld(1:nbndlw+1,0:0,k_topo+1:nzt_rad+1)
             rrtm_lw_tauaer_dum =                                              &
                             rrtm_lw_tauaer(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1)

             CALL rrtmg_lw( 1,                                                 &
                            nzt_rad-k_topo,                                    &
                            rrtm_icld,                                         &
                            rrtm_idrv,                                         &
                            rrtm_play(:,k_topo+1:),                   &
                            rrtm_plev(:,k_topo+1:),                   &
                            rrtm_tlay(:,k_topo+1:),                   &
                            rrtm_tlev(:,k_topo+1:),                   &
                            rrtm_tsfc,                                         &
                            rrtm_h2ovmr(:,k_topo+1:),                 &
                            rrtm_o3vmr(:,k_topo+1:),                  &
                            rrtm_co2vmr(:,k_topo+1:),                 &
                            rrtm_ch4vmr(:,k_topo+1:),                 &
                            rrtm_n2ovmr(:,k_topo+1:),                 &
                            rrtm_o2vmr(:,k_topo+1:),                  &
                            rrtm_cfc11vmr(:,k_topo+1:),               &
                            rrtm_cfc12vmr(:,k_topo+1:),               &
                            rrtm_cfc22vmr(:,k_topo+1:),               &
                            rrtm_ccl4vmr(:,k_topo+1:),                &
                            rrtm_emis,                                         &
                            rrtm_inflglw,                                      &
                            rrtm_iceflglw,                                     &
                            rrtm_liqflglw,                                     &
                            rrtm_cldfr(:,k_topo+1:),                  &
                            rrtm_lw_taucld_dum,                                &
                            rrtm_cicewp(:,k_topo+1:),                 &
                            rrtm_cliqwp(:,k_topo+1:),                 &
                            rrtm_reice(:,k_topo+1:),                  &
                            rrtm_reliq(:,k_topo+1:),                  &
                            rrtm_lw_tauaer_dum,                                &
                            rrtm_lwuflx(:,k_topo:),                   &
                            rrtm_lwdflx(:,k_topo:),                   &
                            rrtm_lwhr(:,k_topo+1:),                   &
                            rrtm_lwuflxc(:,k_topo:),                  &
                            rrtm_lwdflxc(:,k_topo:),                  &
                            rrtm_lwhrc(:,k_topo+1:),                  &
                            rrtm_lwuflx_dt(:,k_topo:),                &
                            rrtm_lwuflxc_dt(:,k_topo:) )

             DEALLOCATE ( rrtm_lw_taucld_dum )
             DEALLOCATE ( rrtm_lw_tauaer_dum )
!
!--          Save fluxes
             DO  i = nxl, nxr
                DO  j = nys, nyn
                   k_topo_l = topo_top_ind(j,i,0)
                   DO k = k_topo_l, nzt+1
                      rad_lw_in(k,j,i)  = rrtm_lwdflx(0,k-k_topo_l+k_topo)
                      rad_lw_out(k,j,i) = rrtm_lwuflx(0,k-k_topo_l+k_topo)
                   ENDDO
                ENDDO
             ENDDO
             rad_lw_in_diff(:,:) = rrtm_lwdflx(0,k_topo)
!
!--          Save heating rates (convert from K/d to K/h).
!--          Further, even though an aggregated radiation is computed, map
!--          signle-column profiles on top of any topography, in order to
!--          obtain correct near surface radiation heating/cooling rates.
             DO  i = nxl, nxr
                DO  j = nys, nyn
                   k_topo_l = topo_top_ind(j,i,0)
                   DO k = k_topo_l+1, nzt+1
                      rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k-k_topo_l+k_topo)  * d_hours_day
                      rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k-k_topo_l+k_topo) * d_hours_day
                   ENDDO
                ENDDO
             ENDDO

          ENDIF

          IF ( sw_radiation .AND. sun_up )  THEN
!
!--          Due to technical reasons, copy optical depths and other
!--          to dummy arguments which are allocated on the exact size as the
!--          rrtmg_sw is called.
!--          As one dimesion is allocated with zero size, compiler complains
!--          that rank of the array does not match that of the
!--          assumed-shaped arguments in the RRTMG library. In order to
!--          avoid this, write to dummy arguments and give pass the entire
!--          dummy array. Seems to be the only existing work-around.
             ALLOCATE( rrtm_sw_taucld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
             ALLOCATE( rrtm_sw_ssacld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
             ALLOCATE( rrtm_sw_asmcld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
             ALLOCATE( rrtm_sw_fsfcld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
             ALLOCATE( rrtm_sw_tauaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
             ALLOCATE( rrtm_sw_ssaaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
             ALLOCATE( rrtm_sw_asmaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
             ALLOCATE( rrtm_sw_ecaer_dum(0:0,k_topo+1:nzt_rad+1,1:naerec+1)  )

             rrtm_sw_taucld_dum = rrtm_sw_taucld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
             rrtm_sw_ssacld_dum = rrtm_sw_ssacld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
             rrtm_sw_asmcld_dum = rrtm_sw_asmcld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
             rrtm_sw_fsfcld_dum = rrtm_sw_fsfcld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
             rrtm_sw_tauaer_dum = rrtm_sw_tauaer(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1)
             rrtm_sw_ssaaer_dum = rrtm_sw_ssaaer(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1)
             rrtm_sw_asmaer_dum = rrtm_sw_asmaer(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1)
             rrtm_sw_ecaer_dum  = rrtm_sw_ecaer(0:0,k_topo+1:nzt_rad+1,1:naerec+1)

             CALL rrtmg_sw( 1,                                                 &
                            nzt_rad-k_topo,                                    &
                            rrtm_icld,                                         &
                            rrtm_iaer,                                         &
                            rrtm_play(:,k_topo+1:nzt_rad+1),                   &
                            rrtm_plev(:,k_topo+1:nzt_rad+2),                   &
                            rrtm_tlay(:,k_topo+1:nzt_rad+1),                   &
                            rrtm_tlev(:,k_topo+1:nzt_rad+2),                   &
                            rrtm_tsfc,                                         &
                            rrtm_h2ovmr(:,k_topo+1:nzt_rad+1),                 &
                            rrtm_o3vmr(:,k_topo+1:nzt_rad+1),                  &
                            rrtm_co2vmr(:,k_topo+1:nzt_rad+1),                 &
                            rrtm_ch4vmr(:,k_topo+1:nzt_rad+1),                 &
                            rrtm_n2ovmr(:,k_topo+1:nzt_rad+1),                 &
                            rrtm_o2vmr(:,k_topo+1:nzt_rad+1),                  &
                            rrtm_asdir,                                        &
                            rrtm_asdif,                                        &
                            rrtm_aldir,                                        &
                            rrtm_aldif,                                        &
                            zenith,                                            &
                            0.0_wp,                                            &
                            day_of_year,                                       &
                            solar_constant,                                    &
                            rrtm_inflgsw,                                      &
                            rrtm_iceflgsw,                                     &
                            rrtm_liqflgsw,                                     &
                            rrtm_cldfr(:,k_topo+1:nzt_rad+1),                  &
                            rrtm_sw_taucld_dum,                                &
                            rrtm_sw_ssacld_dum,                                &
                            rrtm_sw_asmcld_dum,                                &
                            rrtm_sw_fsfcld_dum,                                &
                            rrtm_cicewp(:,k_topo+1:nzt_rad+1),                 &
                            rrtm_cliqwp(:,k_topo+1:nzt_rad+1),                 &
                            rrtm_reice(:,k_topo+1:nzt_rad+1),                  &
                            rrtm_reliq(:,k_topo+1:nzt_rad+1),                  &
                            rrtm_sw_tauaer_dum,                                &
                            rrtm_sw_ssaaer_dum,                                &
                            rrtm_sw_asmaer_dum,                                &
                            rrtm_sw_ecaer_dum,                                 &
                            rrtm_swuflx(:,k_topo:nzt_rad+1),                   &
                            rrtm_swdflx(:,k_topo:nzt_rad+1),                   &
                            rrtm_swhr(:,k_topo+1:nzt_rad+1),                   &
                            rrtm_swuflxc(:,k_topo:nzt_rad+1),                  &
                            rrtm_swdflxc(:,k_topo:nzt_rad+1),                  &
                            rrtm_swhrc(:,k_topo+1:nzt_rad+1),                  &
                            rrtm_dirdflux(:,k_topo:nzt_rad+1),                 &
                            rrtm_difdflux(:,k_topo:nzt_rad+1) )

             DEALLOCATE( rrtm_sw_taucld_dum )
             DEALLOCATE( rrtm_sw_ssacld_dum )
             DEALLOCATE( rrtm_sw_asmcld_dum )
             DEALLOCATE( rrtm_sw_fsfcld_dum )
             DEALLOCATE( rrtm_sw_tauaer_dum )
             DEALLOCATE( rrtm_sw_ssaaer_dum )
             DEALLOCATE( rrtm_sw_asmaer_dum )
             DEALLOCATE( rrtm_sw_ecaer_dum )

!
!--          Save radiation fluxes for the entire depth of the model domain
             DO  i = nxl, nxr
                DO  j = nys, nyn
                   k_topo_l = topo_top_ind(j,i,0)
                   DO k = k_topo_l, nzt+1
                      rad_sw_in(k,j,i)  = rrtm_swdflx(0,k-k_topo_l+k_topo)
                      rad_sw_out(k,j,i) = rrtm_swuflx(0,k-k_topo_l+k_topo)
                   ENDDO
                ENDDO
             ENDDO
!--          Save direct and diffuse SW radiation at the surface (required by RTM)
             rad_sw_in_dir(:,:) = rrtm_dirdflux(0,k_topo)
             rad_sw_in_diff(:,:) = rrtm_difdflux(0,k_topo)

!
!--          Save heating rates (convert from K/d to K/s)
             DO  i = nxl, nxr
                DO  j = nys, nyn
                   k_topo_l = topo_top_ind(j,i,0)
                   DO k = k_topo_l+1, nzt+1
                      rad_sw_hr(k,j,i)     = rrtm_swhr(0,k-k_topo_l+k_topo)  * d_hours_day
                      rad_sw_cs_hr(k,j,i)  = rrtm_swhrc(0,k-k_topo_l+k_topo) * d_hours_day
                   ENDDO
                ENDDO
             ENDDO
!
!--       Solar radiation is zero during night
          ELSE
             rad_sw_in  = 0.0_wp
             rad_sw_out = 0.0_wp
             rad_sw_in_dir(:,:) = 0.0_wp
             rad_sw_in_diff(:,:) = 0.0_wp
          ENDIF
!
!--    RRTMG is called for each (j,i) grid point separately, starting at the
!--    highest topography level. Here no RTM is used since average_radiation is false.
!--    In fact, this branch is only called for homogeneous flat terrain so that the topography-top
!--    index is actually alwasy zero.
       ELSE
!
!--       Loop over all grid points
          DO i = nxl, nxr
             DO j = nys, nyn

!
!--             Prepare profiles of temperature and H2O volume mixing ratio
                DO  m = surf_lsm_h(0)%start_index(j,i), surf_lsm_h(0)%end_index(j,i)
                   rrtm_tlev(0,nzb+1) = surf_lsm_h(0)%pt_surface(m) * exner(nzb)
                ENDDO
                DO  m = surf_usm_h(0)%start_index(j,i), surf_usm_h(0)%end_index(j,i)
                   rrtm_tlev(0,nzb+1) = surf_usm_h(0)%pt_surface(m) * exner(nzb)
                ENDDO


                IF ( bulk_cloud_model )  THEN
                   DO k = nzb+1, nzt+1
                      rrtm_tlay(0,k) = pt(k,j,i) * exner(k)                    &
                                        + lv_d_cp * ql(k,j,i)
                      rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
                   ENDDO
                ELSEIF ( cloud_droplets )  THEN
                   DO k = nzb+1, nzt+1
                      rrtm_tlay(0,k) = pt(k,j,i) * exner(k)                     &
                                        + lv_d_cp * ql(k,j,i)
                      rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q(k,j,i)
                   ENDDO
                ELSE
                   DO k = nzb+1, nzt+1
                      rrtm_tlay(0,k) = pt(k,j,i) * exner(k)
                   ENDDO

                   IF ( humidity )  THEN
                      DO k = nzb+1, nzt+1
                         rrtm_h2ovmr(0,k) =  mol_mass_air_d_wv * q(k,j,i)
                      ENDDO
                   ELSE
                      rrtm_h2ovmr(0,nzb+1:nzt+1) = 0.0_wp
                   ENDIF
                ENDIF

!
!--             Avoid temperature/humidity jumps at the top of the LES domain by
!--             linear interpolation from nzt+2 to nzt+7
                DO k = nzt+2, nzt+7
                   rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1)                         &
                                 + ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) ) &
                                 / ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) ) &
                                 * ( rrtm_play(0,k)     - rrtm_play(0,nzt+1) )

                   rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1)                     &
                              + ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
                              / ( rrtm_play(0,nzt+8)   - rrtm_play(0,nzt+1)   )&
                              * ( rrtm_play(0,k)       - rrtm_play(0,nzt+1) )

                ENDDO

!--             Linear interpolate to zw grid
                DO k = nzb+2, nzt+8
                   rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) -     &
                                      rrtm_tlay(0,k-1))                        &
                                      / ( rrtm_play(0,k) - rrtm_play(0,k-1) )  &
                                      * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
                ENDDO


!
!--             Calculate liquid water path and cloud fraction for each column.
!--             Note that LWP is required in g/m2 instead of kg/kg m.
                rrtm_cldfr  = 0.0_wp
                rrtm_reliq  = 0.0_wp
                rrtm_cliqwp = 0.0_wp
                rrtm_icld   = 0

                IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
                   DO k = nzb+1, nzt+1
                      rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *              &
                                          (rrtm_plev(0,k) - rrtm_plev(0,k+1))  &
                                          * 100.0_wp / g

                      IF ( rrtm_cliqwp(0,k) > 0.0_wp )  THEN
                         rrtm_cldfr(0,k) = 1.0_wp
                         IF ( rrtm_icld == 0 )  rrtm_icld = 1

!
!--                      Calculate cloud droplet effective radius
                         IF ( bulk_cloud_model )  THEN
!
!--                         Calculete effective droplet radius. In case of using
!--                         cloud_scheme = 'morrison' and