Changeset 4624 for palm


Ignore:
Timestamp:
Jul 24, 2020 9:53:17 AM (4 months ago)
Author:
raasch
Message:

file re-formatted to follow the PALM coding standard

File:
1 edited

Legend:

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

    r4587 r4624  
    11!> @file radiation_model_mod.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
    16 !
    17 ! Copyright 2015-2020 Institute of Computer Science of the
    18 !                     Czech Academy of Sciences, Prague
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
     15!
     16! Copyright 2015-2020 Institute of Computer Science of the Czech Academy of Sciences, Prague
    1917! Copyright 2015-2020 Czech Technical University in Prague
    2018! Copyright 1997-2020 Leibniz Universitaet Hannover
    21 !------------------------------------------------------------------------------!
     19!--------------------------------------------------------------------------------------------------!
     20!
    2221!
    2322! Current revisions:
    24 ! ------------------
     23! -----------------
    2524!
    2625!
     
    2827! -----------------
    2928! $Id$
     29! File re-formatted to follow the PALM coding standard
     30!
     31! 4587 2020-07-06 08:53:45Z pavelkrc
    3032! RTM version 3.1 (see also previous commits):
    31 ! - rotation_angle supported
    32 ! - plant canopy box count minimized
    33 ! - multiple enhancements and bugfixes
    34 ! 
     33! - Rotation_angle supported
     34! - Plant canopy box count minimized
     35! - Multiple enhancements and bugfixes
     36!
    3537! 4584 2020-06-29 13:16:14Z pavelkrc
    3638! Consider only boxes with LAD>0 as plant canopy (credit: S. Schubert)
    37 ! 
     39!
    3840! 4576 2020-06-24 17:58:55Z pavelkrc
    3941! Allow the use of rotation_angle in RTM
    40 ! 
     42!
    4143! 4574 2020-06-24 16:33:32Z pavelkrc
    4244! - Restructure code in radiation_check_data_output
    4345! - Move calculation of MPI global array offsets to a subroutine
    44 ! 
     46!
    4547! 4571 2020-06-24 08:59:06Z sebschub
    4648! Bugfix in vertical lad_s coordinate
    47 ! 
     49!
    4850! 4558 2020-06-10 16:27:30Z moh.hefny
    49 ! Bugfix: - reset RTM output average values after each averaging timestep to zero
    50 !         - correct calculation of rtm_rad_net_av
    51 ! 
     51! Bugfix: - Reset RTM output average values after each averaging timestep to zero
     52!         - Correct calculation of rtm_rad_net_av
     53!
    5254! 4555 2020-06-05 21:52:00Z moh.hefny
    5355! Bugfix in averaging PC and MRT related quantities
    54 ! 
     56!
    5557! 4552 2020-06-02 20:33:29Z moh.hefny
    5658! Bugfix in IF statement in the emissivity coupling parameter for radiation-RTM
    57 ! 
     59!
    5860! 4535 2020-05-15 12:07:23Z raasch
    59 ! bugfix for restart data format query
    60 ! 
     61! Bugfix for restart data format query
     62!
    6163! 4531 2020-05-13 09:52:22Z moh.hefny
    6264! Bugfix in gather flux pabs_pc_lwdif in non_parallel case
    6365!
    6466! 4529 2020-05-12 09:14:57Z moh.hefny
    65 ! - added the following new features to the coupling of RTM-radiation model:
    66 !   1) considering the vegetation interaction with LW in the coupling
    67 !   2) considering PC emissivity in calculating the effective emissivity
    68 !   3) new algorithm for claculating the coupling parameters so that each term
    69 !      is calculated within its original line and not at the end.
    70 ! - minor formatting and comments changes
    71 ! 
     67! - Added the following new features to the coupling of RTM-radiation model:
     68!   1) Considering the vegetation interaction with LW in the coupling
     69!   2) Considering PC emissivity in calculating the effective emissivity
     70!   3) New algorithm for claculating the coupling parameters so that each term is calculated within
     71!      its original line and not at the end.
     72! - Minor formatting and comments changes
     73!
    7274! 4517 2020-05-03 14:29:30Z raasch
    73 ! added restart with MPI-IO for reading local arrays
    74 ! 
     75! Added restart with MPI-IO for reading local arrays
     76!
    7577! 4495 2020-04-13 20:11:20Z raasch
    76 ! restart data handling with MPI-IO added
    77 ! 
     78! Restart data handling with MPI-IO added
     79!
    7880! 4493 2020-04-10 09:49:43Z pavelkrc
    7981! Avoid unstable direct normal radiation near horizon
    80 ! 
     82!
    8183! 4481 2020-03-31 18:55:54Z maronga
    82 ! use statement for exchange horiz added
    83 ! 
     84! Use statement for exchange horiz added
     85!
    8486! 4452 2020-03-10 20:15:32Z suehring
    8587! Bugfix in calc_albedo
    8688!
    8789! 4442 2020-03-04 19:21:13Z suehring
    88 ! - Change order of dimension in surface arrays %frac, %emissivity and %albedo
    89 !   to allow for better vectorization in the radiation interactions.
     90! - Change order of dimension in surface arrays %frac, %emissivity and %albedo to allow for better
     91!   vectorization in the radiation interactions.
    9092! - Minor formatting issues
    9193!
    9294! 4441 2020-03-04 19:20:35Z suehring
    93 ! bugfixes: cpp-directives for serial mode moved, small changes to get serial mode compiled
     95! Bugfixes: cpp-directives for serial mode moved, small changes to get serial mode compiled
    9496!
    9597! 4400 2020-02-10 20:32:41Z suehring
     
    102104!
    103105! 4360 2020-01-07 11:25:50Z suehring
    104 ! Renamed pc_heating_rate, pc_transpiration_rate, pc_transpiration_rate to
    105 ! pcm_heating_rate, pcm_latent_rate, pcm_transpiration_rate
     106! Renamed pc_heating_rate, pc_transpiration_rate, pc_transpiration_rate to pcm_heating_rate,
     107! pcm_latent_rate, pcm_transpiration_rate
    106108!
    107109! 4340 2019-12-16 08:17:03Z Giersch
    108 ! Albedo indices for building_surface_pars are now declared as parameters to
    109 ! prevent an error if the gfortran compiler with -Werror=unused-value is used
     110! Albedo indices for building_surface_pars are now declared as parameters to prevent an error if the
     111! gfortran compiler with -Werror=unused-value is used
    110112!
    111113! 4291 2019-11-11 12:36:54Z moh.hefny
    112 ! Enabled RTM in case of biometeorology even if there is no vertical
    113 ! surfaces or 3D vegetation in the domain
     114! Enabled RTM in case of biometeorology even if there is no vertical surfaces or 3D vegetation in
     115! the domain
    114116!
    115117! 4286 2019-10-30 16:01:14Z resler
     
    128130!
    129131! 4227 2019-09-10 18:04:34Z gronemeier
    130 ! implement new palm_date_time_mod
     132! Implement new palm_date_time_mod
    131133!
    132134! 4226 2019-09-10 17:03:24Z suehring
     
    137139! - Revise steering of splitting diffuse and direct radiation
    138140! - Bugfixes in checks
    139 ! - Optimize mapping of radiation components onto 2D arrays, avoid unnecessary
    140 !   operations
     141! - Optimize mapping of radiation components onto 2D arrays, avoid unnecessary operations
    141142!
    142143! 4208 2019-09-02 09:01:07Z suehring
    143 ! Bugfix in accessing albedo_pars in the clear-sky branch
    144 ! (merge from branch resler)
     144! Bugfix in accessing albedo_pars in the clear-sky branch (merge from branch resler)
    145145!
    146146! 4198 2019-08-29 15:17:48Z gronemeier
     
    151151!
    152152! 4190 2019-08-27 15:42:37Z suehring
    153 ! Implement external radiation forcing also for level-of-detail = 2
    154 ! (horizontally 2D radiation)
     153! Implement external radiation forcing also for level-of-detail = 2 (horizontally 2D radiation)
    155154!
    156155! 4188 2019-08-26 14:15:47Z suehring
     
    158157!
    159158! 4187 2019-08-26 12:43:15Z suehring
    160 ! - Take external radiation from root domain dynamic input if not provided for
    161 !   each nested domain
     159! - Take external radiation from root domain dynamic input if not provided for each nested domain
    162160! - Combine MPI_ALLREDUCE calls to reduce mpi overhead
    163161!
     
    187185!
    188186! 4089 2019-07-11 14:30:27Z suehring
    189 ! - Correct level 2 initialization of spectral albedos in rrtmg branch, long- and
    190 !   shortwave albedos were mixed-up.
    191 ! - Change order of albedo_pars so that it is now consistent with the defined
    192 !   order of albedo_pars in PIDS
     187! - Correct level 2 initialization of spectral albedos in rrtmg branch, long- and shortwave albedos
     188!   were mixed-up.
     189! - Change order of albedo_pars so that it is now consistent with the defined order of albedo_pars
     190!   in PIDS
    193191!
    194192! 4069 2019-07-01 14:05:51Z Giersch
    195 ! Masked output running index mid has been introduced as a local variable to
    196 ! avoid runtime error (Loop variable has been modified) in time_integration
     193! Masked output running index mid has been introduced as a local variable to avoid runtime error
     194! (Loop variable has been modified) in time_integration
    197195!
    198196! 4067 2019-07-01 13:29:25Z suehring
     
    203201!
    204202! 4008 2019-05-30 09:50:11Z moh.hefny
    205 ! Bugfix in check variable when a variable's string is less than 3
    206 ! characters is processed. All variables now are checked if they
    207 ! belong to radiation
     203! Bugfix in check variable when a variable's string is less than 3 characters is processed. All
     204! variables now are checked if they belong to radiation
    208205!
    209206! 3992 2019-05-22 16:49:38Z suehring
    210 ! Bugfix in rrtmg radiation branch in a nested run when the lowest prognistic
    211 ! grid points in a child domain are all inside topography
     207! Bugfix in rrtmg radiation branch in a nested run when the lowest prognistic grid points in a child
     208! domain are all inside topography
    212209!
    213210! 3987 2019-05-22 09:52:13Z kanani
     
    221218!
    222219! 3885 2019-04-11 11:29:34Z kanani
    223 ! Changes related to global restructuring of location messages and introduction
    224 ! of additional debug messages
     220! Changes related to global restructuring of location messages and introduction of additional debug
     221! messages
    225222!
    226223! 3881 2019-04-10 09:31:22Z suehring
    227 ! Output of albedo and emissivity moved from USM, bugfixes in initialization
    228 ! of albedo
     224! Output of albedo and emissivity moved from USM, bugfixes in initialization of albedo
    229225!
    230226! 3861 2019-04-04 06:27:41Z maronga
     
    238234!
    239235! 3846 2019-04-01 13:55:30Z suehring
    240 ! unused variable removed
     236! Unused variable removed
    241237!
    242238! 3814 2019-03-26 08:40:31Z pavelkrc
     
    246242!
    247243! 3771 2019-02-28 12:19:33Z raasch
    248 ! rrtmg preprocessor for directives moved/added, save attribute added to temporary
    249 ! pointers to avoid compiler warnings about outlived pointer targets,
    250 ! statement added to avoid compiler warning about unused variable
     244! rrtmg preprocessor for directives moved/added, save attribute added to temporary pointers to avoid
     245! compiler warnings about outlived pointer targets, statement added to avoid compiler warning about
     246! unused variable
    251247!
    252248! 3769 2019-02-28 10:16:49Z moh.hefny
    253 ! removed unused variables and subroutine radiation_radflux_gridbox
     249! Removed unused variables and subroutine radiation_radflux_gridbox
    254250!
    255251! 3767 2019-02-27 08:18:02Z raasch
    256 ! unused variable for file index removed from rrd-subroutines parameter list
     252! Unused variable for file index removed from rrd-subroutines parameter list
    257253!
    258254! 3760 2019-02-21 18:47:35Z moh.hefny
    259 ! Bugfix: initialized simulated_time before calculating solar position
    260 ! to enable restart option with reading in SVF from file(s).
     255! Bugfix: initialized simulated_time before calculating solar position to enable restart option with
     256! reading in SVF from file(s).
    261257!
    262258! 3754 2019-02-19 17:02:26Z kanani
    263259! (resler, pavelkrc)
    264 ! Bugfixes: add further required MRT factors to read/write_svf,
    265 ! fix for aggregating view factors to eliminate local noise in reflected
    266 ! irradiance at mutually close surfaces (corners, presence of trees) in the
    267 ! angular discretization scheme.
     260! Bugfixes: add further required MRT factors to read/write_svf, fix for aggregating view factors to
     261! eliminate local noise in reflected irradiance at mutually close surfaces (corners, presence of
     262! trees) in the angular discretization scheme.
    268263!
    269264! 3752 2019-02-19 09:37:22Z resler
    270 ! added read/write number of MRT factors to the respective routines
     265! Added read/write number of MRT factors to the respective routines
    271266!
    272267! 3705 2019-01-29 19:56:39Z suehring
     
    280275!
    281276! 3655 2019-01-07 16:51:22Z knoop
    282 ! nopointer option removed
     277! Nopointer option removed
    283278!
    284279! 1496 2014-12-02 17:25:50Z maronga
     
    286281!
    287282!
     283!--------------------------------------------------------------------------------------------------!
    288284! Description:
    289285! ------------
    290286!> Radiation models and interfaces:
    291 !> constant, simple and RRTMG models, interface to external radiation model
     287!> Constant, simple and RRTMG models, interface to external radiation model
    292288!> Radiative Transfer Model (RTM) version 3.0 for modelling of radiation
    293 !> interactions within urban canopy or other surface layer in complex terrain
     289!> Interactions within urban canopy or other surface layer in complex terrain
    294290!> Integrations of RTM with other PALM-4U modules:
    295 !> integration with RRTMG, USM, LSM, PCM, BIO modules
     291!> Integration with RRTMG, USM, LSM, PCM, BIO modules
    296292!>
    297 !> @todo move variable definitions used in radiation_init only to the subroutine
    298 !>       as they are no longer required after initialization.
     293!> @todo Move variable definitions used in radiation_init only to the subroutine as they are no
     294!>       longer required after initialization.
    299295!> @todo Output of full column vertical profiles used in RRTMG
    300296!> @todo Output of other rrtm arrays (such as volume mixing ratios)
    301297!> @todo Optimize radiation_tendency routines
    302298!>
    303 !> @note Many variables have a leading dummy dimension (0:0) in order to
    304 !>       match the assume-size shape expected by the RRTMG model.
    305 !------------------------------------------------------------------------------!
     299!> @note Many variables have a leading dummy dimension (0:0) in order to match the assume-size shape
     300!>       expected by the RRTMG model.
     301!--------------------------------------------------------------------------------------------------!
    306302 MODULE radiation_model_mod
    307303
    308     USE arrays_3d,                                                             &
    309         ONLY:  dzw, hyp, nc, pt, p, q, ql, u, v, w, zu, zw, exner, d_exner
    310 
    311     USE basic_constants_and_equations_mod,                                     &
    312         ONLY:  c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, sigma_sb, &
    313                barometric_formula
    314 
    315     USE calc_mean_profile_mod,                                                 &
     304    USE arrays_3d,                                                                                 &
     305        ONLY:  dzw,                                                                                &
     306               d_exner,                                                                            &
     307               exner,                                                                              &
     308               hyp,                                                                                &
     309               nc,                                                                                 &
     310               pt,                                                                                 &
     311               p,                                                                                  &
     312               q,                                                                                  &
     313               ql,                                                                                 &
     314               u,                                                                                  &
     315               v,                                                                                  &
     316               w,                                                                                  &
     317               zu,                                                                                 &
     318               zw
     319
     320
     321
     322    USE basic_constants_and_equations_mod,                                                         &
     323        ONLY:  barometric_formula,                                                                 &
     324               c_p,                                                                                &
     325               g,                                                                                  &
     326               lv_d_cp,                                                                            &
     327               l_v,                                                                                &
     328               pi,                                                                                 &
     329               r_d,                                                                                &
     330               rho_l,                                                                              &
     331               sigma_sb,                                                                           &
     332               solar_constant
     333
     334
     335
     336    USE calc_mean_profile_mod,                                                                     &
    316337        ONLY:  calc_mean_profile
    317338
    318     USE control_parameters,                                                    &
    319         ONLY:  biometeorology, cloud_droplets, coupling_char,                  &
    320                debug_output, debug_output_timestep, debug_string,              &
    321                dt_3d,                                                          &
    322                dz, dt_spinup, end_time,                                        &
    323                humidity,                                                       &
    324                initializing_actions, io_blocks, io_group,                      &
    325                land_surface, large_scale_forcing,                              &
    326                latitude, longitude, lsf_surf,                                  &
    327                message_string, plant_canopy, pt_surface,                       &
    328                rho_surface, simulated_time, spinup_time, surface_pressure,     &
    329                read_svf, restart_data_format_output, write_svf,                &
    330                time_since_reference_point, urban_surface, varnamelength
    331 
    332     USE cpulog,                                                                &
    333         ONLY:  cpu_log, log_point, log_point_s
    334 
    335     USE grid_variables,                                                        &
    336          ONLY:  ddx, ddy, dx, dy
    337 
    338     USE indices,                                                               &
    339         ONLY:  nnx, nny, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,   &
    340                nzb, nzt, topo_top_ind
     339    USE control_parameters,                                                                        &
     340        ONLY:  biometeorology,                                                                     &
     341               cloud_droplets,                                                                     &
     342               coupling_char,                                                                      &
     343               debug_output,                                                                       &
     344               debug_output_timestep,                                                              &
     345               debug_string,                                                                       &
     346               dt_3d,                                                                              &
     347               dz,                                                                                 &
     348               dt_spinup,                                                                          &
     349               end_time,                                                                           &
     350               humidity,                                                                           &
     351               initializing_actions,                                                               &
     352               io_blocks,                                                                          &
     353               io_group,                                                                           &
     354               land_surface,                                                                       &
     355               large_scale_forcing,                                                                &
     356               latitude,                                                                           &
     357               longitude,                                                                          &
     358               lsf_surf,                                                                           &
     359               message_string,                                                                     &
     360               plant_canopy,                                                                       &
     361               pt_surface,                                                                         &
     362               read_svf,                                                                           &
     363               restart_data_format_output,                                                         &
     364               rho_surface,                                                                        &
     365               simulated_time,                                                                     &
     366               spinup_time,                                                                        &
     367               surface_pressure,                                                                   &
     368               time_since_reference_point,                                                         &
     369               urban_surface,                                                                      &
     370               varnamelength,                                                                      &
     371               write_svf
     372
     373    USE cpulog,                                                                                    &
     374        ONLY:  cpu_log,                                                                            &
     375               log_point,                                                                          &
     376               log_point_s
     377
     378    USE grid_variables,                                                                            &
     379         ONLY:  ddx,                                                                               &
     380                ddy,                                                                               &
     381                dx,                                                                                &
     382                dy
     383
     384    USE indices,                                                                                   &
     385        ONLY:  nnx,                                                                                &
     386               nny,                                                                                &
     387               nx,                                                                                 &
     388               nxl,                                                                                &
     389               nxlg,                                                                               &
     390               nxr,                                                                                &
     391               nxrg,                                                                               &
     392               ny,                                                                                 &
     393               nyn,                                                                                &
     394               nyng,                                                                               &
     395               nys,                                                                                &
     396               nysg,                                                                               &
     397               nzb,                                                                                &
     398               nzt,                                                                                &
     399               topo_top_ind
    341400
    342401    USE, INTRINSIC :: iso_c_binding
     
    344403    USE kinds
    345404
    346     USE bulk_cloud_model_mod,                                                  &
    347         ONLY:  bulk_cloud_model, microphysics_morrison, na_init, nc_const, sigma_gc
     405    USE bulk_cloud_model_mod,                                                                      &
     406        ONLY:  bulk_cloud_model,                                                                   &
     407               microphysics_morrison,                                                              &
     408               na_init,                                                                            &
     409               nc_const,                                                                           &
     410               sigma_gc
    348411
    349412#if defined ( __netcdf )
     
    351414#endif
    352415
    353     USE netcdf_data_input_mod,                                                 &
    354         ONLY:  albedo_type_f,                                                  &
    355                albedo_pars_f,                                                  &
    356                building_type_f,                                                &
    357                building_surface_pars_f,                                        &
    358                pavement_type_f,                                                &
    359                vegetation_type_f,                                              &
    360                water_type_f,                                                   &
    361                char_fill,                                                      &
    362                char_lod,                                                       &
    363                check_existence,                                                &
    364                close_input_file,                                               &
    365                get_attribute,                                                  &
    366                get_dimension_length,                                           &
    367                get_variable,                                                   &
    368                inquire_num_variables,                                          &
    369                inquire_variable_names,                                         &
    370                input_file_dynamic,                                             &
    371                input_pids_dynamic,                                             &
    372                num_var_pids,                                                   &
    373                pids_id,                                                        &
    374                open_read_file,                                                 &
    375                real_1d_3d,                                                     &
    376                vars_pids
    377 
    378     USE palm_date_time_mod,                                                    &
    379         ONLY:  date_time_str_len, get_date_time,                               &
    380                hours_per_day, seconds_per_hour
    381 
    382     USE plant_canopy_model_mod,                                                &
    383         ONLY:  lad_s,                                                          &
    384                pcm_heating_rate,                                               &
    385                pcm_transpiration_rate,                                         &
    386                pcm_latent_rate,                                                &
    387                plant_canopy_transpiration,                                     &
    388                pcm_calc_transpiration_rate
     416    USE netcdf_data_input_mod,                                                                     &
     417        ONLY:  albedo_type_f,                                                                      &
     418               albedo_pars_f,                                                                      &
     419               building_type_f,                                                                    &
     420               building_surface_pars_f,                                                            &
     421               char_fill,                                                                          &
     422               char_lod,                                                                           &
     423               check_existence,                                                                    &
     424               close_input_file,                                                                   &
     425               get_attribute,                                                                      &
     426               get_dimension_length,                                                               &
     427               get_variable,                                                                       &
     428               inquire_num_variables,                                                              &
     429               inquire_variable_names,                                                             &
     430               input_file_dynamic,                                                                 &
     431               input_pids_dynamic,                                                                 &
     432               num_var_pids,                                                                       &
     433               open_read_file,                                                                     &
     434               pavement_type_f,                                                                    &
     435               pids_id,                                                                            &
     436               real_1d_3d,                                                                         &
     437               vars_pids,                                                                          &
     438               vegetation_type_f,                                                                  &
     439               water_type_f
     440
     441
     442
     443    USE palm_date_time_mod,                                                                        &
     444        ONLY:  date_time_str_len,                                                                  &
     445               get_date_time,                                                                      &
     446               hours_per_day,                                                                      &
     447               seconds_per_hour
     448
     449    USE plant_canopy_model_mod,                                                                    &
     450        ONLY:  lad_s,                                                                              &
     451               pcm_calc_transpiration_rate,                                                        &
     452               pcm_heating_rate,                                                                   &
     453               pcm_transpiration_rate,                                                             &
     454               pcm_latent_rate,                                                                    &
     455               plant_canopy_transpiration
     456
    389457
    390458    USE pegrid
    391459
    392460#if defined ( __rrtmg )
    393     USE parrrsw,                                                               &
    394         ONLY:  naerec, nbndsw
    395 
    396     USE parrrtm,                                                               &
     461    USE parrrsw,                                                                                   &
     462        ONLY:  naerec,                                                                             &
     463               nbndsw
     464
     465    USE parrrtm,                                                                                   &
    397466        ONLY:  nbndlw
    398467
    399     USE rrtmg_lw_init,                                                         &
     468    USE rrtmg_lw_init,                                                                             &
    400469        ONLY:  rrtmg_lw_ini
    401470
    402     USE rrtmg_sw_init,                                                         &
     471    USE rrtmg_sw_init,                                                                             &
    403472        ONLY:  rrtmg_sw_ini
    404473
    405     USE rrtmg_lw_rad,                                                          &
     474    USE rrtmg_lw_rad,                                                                              &
    406475        ONLY:  rrtmg_lw
    407476
    408     USE rrtmg_sw_rad,                                                          &
     477    USE rrtmg_sw_rad,                                                                              &
    409478        ONLY:  rrtmg_sw
    410479#endif
    411480    USE restart_data_mpi_io_mod,                                                                   &
    412         ONLY:  rd_mpi_io_check_array, rrd_mpi_io, wrd_mpi_io
    413 
    414     USE statistics,                                                            &
     481        ONLY:  rd_mpi_io_check_array,                                                              &
     482               rrd_mpi_io,                                                                         &
     483               wrd_mpi_io
     484
     485    USE statistics,                                                                                &
    415486        ONLY:  hom
    416487
    417     USE surface_mod,                                                           &
    418         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win,                       &
    419                surf_lsm_h, surf_lsm_v, surf_type, surf_usm_h, surf_usm_v,      &
     488    USE surface_mod,                                                                               &
     489        ONLY:  ind_pav_green,                                                                      &
     490               ind_veg_wall,                                                                       &
     491               ind_wat_win,                                                                        &
     492               surf_lsm_h,                                                                         &
     493               surf_lsm_v,                                                                         &
     494               surf_type,                                                                          &
     495               surf_usm_h,                                                                         &
     496               surf_usm_v,                                                                         &
    420497               vertical_surfaces_exist
    421498
    422499    IMPLICIT NONE
    423500
    424     CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
     501    CHARACTER(10) ::  radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
    425502
    426503!
    427504!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
    428     CHARACTER(37), DIMENSION(0:33), PARAMETER :: albedo_type_name = (/      &
    429                                    'user defined                         ', & !  0
    430                                    'ocean                                ', & !  1
    431                                    'mixed farming, tall grassland        ', & !  2
    432                                    'tall/medium grassland                ', & !  3
    433                                    'evergreen shrubland                  ', & !  4
    434                                    'short grassland/meadow/shrubland     ', & !  5
    435                                    'evergreen needleleaf forest          ', & !  6
    436                                    'mixed deciduous evergreen forest     ', & !  7
    437                                    'deciduous forest                     ', & !  8
    438                                    'tropical evergreen broadleaved forest', & !  9
    439                                    'medium/tall grassland/woodland       ', & ! 10
    440                                    'desert, sandy                        ', & ! 11
    441                                    'desert, rocky                        ', & ! 12
    442                                    'tundra                               ', & ! 13
    443                                    'land ice                             ', & ! 14
    444                                    'sea ice                              ', & ! 15
    445                                    'snow                                 ', & ! 16
    446                                    'bare soil                            ', & ! 17
    447                                    'asphalt/concrete mix                 ', & ! 18
    448                                    'asphalt (asphalt concrete)           ', & ! 19
    449                                    'concrete (Portland concrete)         ', & ! 20
    450                                    'sett                                 ', & ! 21
    451                                    'paving stones                        ', & ! 22
    452                                    'cobblestone                          ', & ! 23
    453                                    'metal                                ', & ! 24
    454                                    'wood                                 ', & ! 25
    455                                    'gravel                               ', & ! 26
    456                                    'fine gravel                          ', & ! 27
    457                                    'pebblestone                          ', & ! 28
    458                                    'woodchips                            ', & ! 29
    459                                    'tartan (sports)                      ', & ! 30
    460                                    'artifical turf (sports)              ', & ! 31
    461                                    'clay (sports)                        ', & ! 32
    462                                    'building (dummy)                     '  & ! 33
     505    CHARACTER(37), DIMENSION(0:33), PARAMETER :: albedo_type_name = (/                       &
     506                                   'user defined                         ',                  &  !<  0
     507                                   'ocean                                ',                  &  !<  1
     508                                   'mixed farming, tall grassland        ',                  &  !<  2
     509                                   'tall/medium grassland                ',                  &  !<  3
     510                                   'evergreen shrubland                  ',                  &  !<  4
     511                                   'short grassland/meadow/shrubland     ',                  &  !<  5
     512                                   'evergreen needleleaf forest          ',                  &  !<  6
     513                                   'mixed deciduous evergreen forest     ',                  &  !<  7
     514                                   'deciduous forest                     ',                  &  !<  8
     515                                   'tropical evergreen broadleaved forest',                  &  !<  9
     516                                   'medium/tall grassland/woodland       ',                  &  !< 10
     517                                   'desert, sandy                        ',                  &  !< 11
     518                                   'desert, rocky                        ',                  &  !< 12
     519                                   'tundra                               ',                  &  !< 13
     520                                   'land ice                             ',                  &  !< 14
     521                                   'sea ice                              ',                  &  !< 15
     522                                   'snow                                 ',                  &  !< 16
     523                                   'bare soil                            ',                  &  !< 17
     524                                   'asphalt/concrete mix                 ',                  &  !< 18
     525                                   'asphalt (asphalt concrete)           ',                  &  !< 19
     526                                   'concrete (Portland concrete)         ',                  &  !< 20
     527                                   'sett                                 ',                  &  !< 21
     528                                   'paving stones                        ',                  &  !< 22
     529                                   'cobblestone                          ',                  &  !< 23
     530                                   'metal                                ',                  &  !< 24
     531                                   'wood                                 ',                  &  !< 25
     532                                   'gravel                               ',                  &  !< 26
     533                                   'fine gravel                          ',                  &  !< 27
     534                                   'pebblestone                          ',                  &  !< 28
     535                                   'woodchips                            ',                  &  !< 29
     536                                   'tartan (sports)                      ',                  &  !< 30
     537                                   'artifical turf (sports)              ',                  &  !< 31
     538                                   'clay (sports)                        ',                  &  !< 32
     539                                   'building (dummy)                     '                   &  !< 33
    463540                                                         /)
    464541!
    465 !-- Indices of radiation-related input attributes in building_surface_pars
    466 !-- (other are in urban_surface_mod)
    467     INTEGER(iwp), PARAMETER ::  ind_s_alb_b_wall                = 19 !< index for Broadband albedo of wall fraction
    468     INTEGER(iwp), PARAMETER ::  ind_s_alb_l_wall                = 20 !< index for Longwave albedo of wall fraction
    469     INTEGER(iwp), PARAMETER ::  ind_s_alb_s_wall                = 21 !< index for Shortwave albedo of wall fraction
    470     INTEGER(iwp), PARAMETER ::  ind_s_alb_b_win                 = 22 !< index for Broadband albedo of window fraction
    471     INTEGER(iwp), PARAMETER ::  ind_s_alb_l_win                 = 23 !< index for Longwave albedo of window fraction
    472     INTEGER(iwp), PARAMETER ::  ind_s_alb_s_win                 = 24 !< index for Shortwave albedo of window fraction
    473     INTEGER(iwp), PARAMETER ::  ind_s_alb_b_green               = 24 !< index for Broadband albedo of green fraction
    474     INTEGER(iwp), PARAMETER ::  ind_s_alb_l_green               = 25 !< index for Longwave albedo of green fraction
    475     INTEGER(iwp), PARAMETER ::  ind_s_alb_s_green               = 26 !< index for Shortwave albedo of green fraction
    476 
    477     INTEGER(iwp) :: albedo_type  = 9999999_iwp, &     !< Albedo surface type
    478                     dots_rad     = 0_iwp              !< starting index for timeseries output
    479 
    480     LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
    481                 constant_albedo = .FALSE.,            & !< flag parameter indicating whether the albedo may change depending on zenith
    482                 force_radiation_call = .FALSE.,       & !< flag parameter for unscheduled radiation calls
    483                 lw_radiation = .TRUE.,                & !< flag parameter indicating whether longwave radiation shall be calculated
    484                 radiation = .FALSE.,                  & !< flag parameter indicating whether the radiation model is used
    485                 sun_up    = .TRUE.,                   & !< flag parameter indicating whether the sun is up or down
    486                 sw_radiation = .TRUE.,                & !< flag parameter indicating whether shortwave radiation shall be calculated
    487                 sun_direction = .FALSE.,              & !< flag parameter indicating whether solar direction shall be calculated
    488                 average_radiation = .FALSE.,          & !< flag to set the calculation of radiation averaging for the domain
    489                 radiation_interactions = .FALSE.,     & !< flag to activiate RTM (TRUE only if vertical urban/land surface and trees exist)
    490                 surface_reflections = .TRUE.,         & !< flag to switch the calculation of radiation interaction between surfaces.
    491                                                         !< When it switched off, only the effect of buildings and trees shadow
    492                                                         !< will be considered. However fewer SVFs are expected.
    493                 radiation_interactions_on = .TRUE.      !< namelist flag to force RTM activiation regardless to vertical urban/land surface and trees
    494 
    495     REAL(wp) :: albedo = 9999999.9_wp,           & !< NAMELIST alpha
    496                 albedo_lw_dif = 9999999.9_wp,    & !< NAMELIST aldif
    497                 albedo_lw_dir = 9999999.9_wp,    & !< NAMELIST aldir
    498                 albedo_sw_dif = 9999999.9_wp,    & !< NAMELIST asdif
    499                 albedo_sw_dir = 9999999.9_wp,    & !< NAMELIST asdir
    500                 decl_1,                          & !< declination coef. 1
    501                 decl_2,                          & !< declination coef. 2
    502                 decl_3,                          & !< declination coef. 3
    503                 dt_radiation = 0.0_wp,           & !< radiation model timestep
    504                 emissivity = 9999999.9_wp,       & !< NAMELIST surface emissivity
    505                 lon = 0.0_wp,                    & !< longitude in radians
    506                 lat = 0.0_wp,                    & !< latitude in radians
    507                 net_radiation = 0.0_wp,          & !< net radiation at surface
    508                 skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
    509                 sky_trans,                       & !< sky transmissivity
    510                 time_radiation = 0.0_wp,         & !< time since last call of radiation code
    511                 trace_fluxes_above = -1.0_wp,    & !< NAMELIST option for debug tracing of large radiative fluxes (W/m2;W/m3)
    512                 min_stable_coszen = 0.0262_wp      !< 1.5 deg above horizon, eliminates most of circumsolar
    513 
    514     INTEGER(iwp) ::  day_of_year   !< day of the current year
    515 
    516     REAL(wp) ::  cos_zenith        !< cosine of solar zenith angle, also z-coordinate of solar unit vector
    517     REAL(wp) ::  d_hours_day       !< 1 / hours-per-day
    518     REAL(wp) ::  d_seconds_hour    !< 1 / seconds-per-hour
    519     REAL(wp) ::  second_of_day     !< second of the current day
    520     REAL(wp) ::  sun_dir_lat       !< y-coordinate of solar unit vector
    521     REAL(wp) ::  sun_dir_lon       !< x-coordinate of solar unit vector
    522 
    523     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_net_av       !< average of net radiation (rad_net) at surface
    524     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_in_xy_av  !< average of incoming longwave radiation at surface
    525     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_out_xy_av !< average of outgoing longwave radiation at surface
    526     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_in_xy_av  !< average of incoming shortwave radiation at surface
    527     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_out_xy_av !< average of outgoing shortwave radiation at surface
    528 
    529     REAL(wp), PARAMETER :: emissivity_atm_clsky = 0.8_wp       !< emissivity of the clear-sky atmosphere
     542!-- Indices of radiation-related input attributes in building_surface_pars (others are in
     543!-- urban_surface_mod)
     544    INTEGER(iwp), PARAMETER ::  ind_s_alb_b_wall  = 19  !< index for Broadband albedo of wall fraction
     545    INTEGER(iwp), PARAMETER ::  ind_s_alb_l_wall  = 20  !< index for Longwave albedo of wall fraction
     546    INTEGER(iwp), PARAMETER ::  ind_s_alb_s_wall  = 21  !< index for Shortwave albedo of wall fraction
     547    INTEGER(iwp), PARAMETER ::  ind_s_alb_b_win   = 22  !< index for Broadband albedo of window fraction
     548    INTEGER(iwp), PARAMETER ::  ind_s_alb_l_win   = 23  !< index for Longwave albedo of window fraction
     549    INTEGER(iwp), PARAMETER ::  ind_s_alb_s_win   = 24  !< index for Shortwave albedo of window fraction
     550    INTEGER(iwp), PARAMETER ::  ind_s_alb_b_green = 24  !< index for Broadband albedo of green fraction
     551    INTEGER(iwp), PARAMETER ::  ind_s_alb_l_green = 25  !< index for Longwave albedo of green fraction
     552    INTEGER(iwp), PARAMETER ::  ind_s_alb_s_green = 26  !< index for Shortwave albedo of green fraction
     553
     554    INTEGER(iwp) ::  albedo_type  = 9999999_iwp, &  !< Albedo surface type
     555                     dots_rad     = 0_iwp           !< starting index for timeseries output
     556    INTEGER(iwp) ::  day_of_year                    !< day of the current year
     557
     558    LOGICAL ::  unscheduled_radiation_calls = .TRUE., &  !< flag parameter indicating whether additional calls of the radiation
     559                                                         !< code are allowed
     560                constant_albedo =            .FALSE., &  !< flag parameter indicating whether the albedo may change depending on
     561                                                         !< zenith
     562                force_radiation_call =       .FALSE., &  !< flag parameter for unscheduled radiation calls
     563                lw_radiation =                .TRUE., &  !< flag parameter indicating whether longwave radiation shall be calculated
     564                radiation =                  .FALSE., &  !< flag parameter indicating whether the radiation model is used
     565                sun_up =                      .TRUE., &  !< flag parameter indicating whether the sun is up or down
     566                sw_radiation =                .TRUE., &  !< flag parameter indicating whether shortwave radiation shall be
     567                                                         !< calculated
     568                sun_direction =              .FALSE., &  !< flag parameter indicating whether solar direction shall be calculated
     569                average_radiation =          .FALSE., &  !< flag to set the calculation of radiation averaging for the domain
     570                radiation_interactions =     .FALSE., &  !< flag to activiate RTM (TRUE only if vertical urban/land surface and
     571                                                         !< trees exist)
     572                surface_reflections =         .TRUE., &  !< flag to switch the calculation of radiation interaction between
     573                                                         !< surfaces. When it switched off, only the effect of buildings and trees
     574                                                         !< shadow will be considered. However fewer SVFs are expected.
     575                radiation_interactions_on = .TRUE.       !< namelist flag to force RTM activiation regardless to vertical urban/
     576                                                         !< land surface and trees
     577
     578    REAL(wp) ::  albedo = 9999999.9_wp,           &  !< NAMELIST alpha
     579                 albedo_lw_dif = 9999999.9_wp,    &  !< NAMELIST aldif
     580                 albedo_lw_dir = 9999999.9_wp,    &  !< NAMELIST aldir
     581                 albedo_sw_dif = 9999999.9_wp,    &  !< NAMELIST asdif
     582                 albedo_sw_dir = 9999999.9_wp,    &  !< NAMELIST asdir
     583                 decl_1,                          &  !< declination coef. 1
     584                 decl_2,                          &  !< declination coef. 2
     585                 decl_3,                          &  !< declination coef. 3
     586                 dt_radiation = 0.0_wp,           &  !< radiation model timestep
     587                 emissivity = 9999999.9_wp,       &  !< NAMELIST surface emissivity
     588                 lon = 0.0_wp,                    &  !< longitude in radians
     589                 lat = 0.0_wp,                    &  !< latitude in radians
     590                 net_radiation = 0.0_wp,          &  !< net radiation at surface
     591                 skip_time_do_radiation = 0.0_wp, &  !< Radiation model is not called before this time
     592                 sky_trans,                       &  !< sky transmissivity
     593                 time_radiation = 0.0_wp,         &  !< time since last call of radiation code
     594                 trace_fluxes_above = -1.0_wp,    &  !< NAMELIST option for debug tracing of large radiative fluxes (W/m2;W/m3)
     595                 min_stable_coszen = 0.0262_wp       !< 1.5 deg above horizon, eliminates most of circumsolar
     596
     597    REAL(wp) ::  cos_zenith      !< cosine of solar zenith angle, also z-coordinate of solar unit vector
     598    REAL(wp) ::  d_hours_day     !< 1 / hours-per-day
     599    REAL(wp) ::  d_seconds_hour  !< 1 / seconds-per-hour
     600    REAL(wp) ::  second_of_day   !< second of the current day
     601    REAL(wp) ::  sun_dir_lat     !< y-coordinate of solar unit vector
     602    REAL(wp) ::  sun_dir_lon     !< x-coordinate of solar unit vector
     603
     604    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_net_av        !< average of net radiation (rad_net) at surface
     605    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_in_xy_av   !< average of incoming longwave radiation at surface
     606    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_out_xy_av  !< average of outgoing longwave radiation at surface
     607    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_in_xy_av   !< average of incoming shortwave radiation at surface
     608    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_out_xy_av  !< average of outgoing shortwave radiation at surface
     609
     610    REAL(wp), PARAMETER ::  emissivity_atm_clsky = 0.8_wp  !< emissivity of the clear-sky atmosphere
    530611!
    531612!-- Land surface albedos for solar zenith angle of 60degree after Briegleb (1992)
    532 !-- (broadband, longwave, shortwave ):   bb,      lw,      sw,
    533     REAL(wp), DIMENSION(0:2,1:33), PARAMETER :: albedo_pars = RESHAPE( (/&
    534                                    0.06_wp, 0.06_wp, 0.06_wp,            & !  1
    535                                    0.19_wp, 0.28_wp, 0.09_wp,            & !  2
    536                                    0.23_wp, 0.33_wp, 0.11_wp,            & !  3
    537                                    0.23_wp, 0.33_wp, 0.11_wp,            & !  4
    538                                    0.25_wp, 0.34_wp, 0.14_wp,            & !  5
    539                                    0.14_wp, 0.22_wp, 0.06_wp,            & !  6
    540                                    0.17_wp, 0.27_wp, 0.06_wp,            & !  7
    541                                    0.19_wp, 0.31_wp, 0.06_wp,            & !  8
    542                                    0.14_wp, 0.22_wp, 0.06_wp,            & !  9
    543                                    0.18_wp, 0.28_wp, 0.06_wp,            & ! 10
    544                                    0.43_wp, 0.51_wp, 0.35_wp,            & ! 11
    545                                    0.32_wp, 0.40_wp, 0.24_wp,            & ! 12
    546                                    0.19_wp, 0.27_wp, 0.10_wp,            & ! 13
    547                                    0.77_wp, 0.65_wp, 0.90_wp,            & ! 14
    548                                    0.77_wp, 0.65_wp, 0.90_wp,            & ! 15
    549                                    0.82_wp, 0.70_wp, 0.95_wp,            & ! 16
    550                                    0.08_wp, 0.08_wp, 0.08_wp,            & ! 17
    551                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 18
    552                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 19
    553                                    0.30_wp, 0.30_wp, 0.30_wp,            & ! 20
    554                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 21
    555                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 22
    556                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 23
    557                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 24
    558                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 25
    559                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 26
    560                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 27
    561                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 28
    562                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 29
    563                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 30
    564                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 31
    565                                    0.17_wp, 0.17_wp, 0.17_wp,            & ! 32
    566                                    0.17_wp, 0.17_wp, 0.17_wp             & ! 33
    567                                  /), (/ 3, 33 /) )
    568 
    569     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
    570                         rad_lw_cs_hr,                  & !< longwave clear sky radiation heating rate (K/s)
    571                         rad_lw_cs_hr_av,               & !< average of rad_lw_cs_hr
    572                         rad_lw_hr,                     & !< longwave radiation heating rate (K/s)
    573                         rad_lw_hr_av,                  & !< average of rad_sw_hr
    574                         rad_lw_in,                     & !< incoming longwave radiation (W/m2)
    575                         rad_lw_in_av,                  & !< average of rad_lw_in
    576                         rad_lw_out,                    & !< outgoing longwave radiation (W/m2)
    577                         rad_lw_out_av,                 & !< average of rad_lw_out
    578                         rad_sw_cs_hr,                  & !< shortwave clear sky radiation heating rate (K/s)
    579                         rad_sw_cs_hr_av,               & !< average of rad_sw_cs_hr
    580                         rad_sw_hr,                     & !< shortwave radiation heating rate (K/s)
    581                         rad_sw_hr_av,                  & !< average of rad_sw_hr
    582                         rad_sw_in,                     & !< incoming shortwave radiation (W/m2)
    583                         rad_sw_in_av,                  & !< average of rad_sw_in
    584                         rad_sw_out,                    & !< outgoing shortwave radiation (W/m2)
    585                         rad_sw_out_av                    !< average of rad_sw_out
     613!-- (broadband, longwave, shortwave ):              bb,      lw,      sw,
     614    REAL(wp), DIMENSION(0:2,1:33), PARAMETER ::  albedo_pars = RESHAPE( (/&
     615                                                 0.06_wp, 0.06_wp, 0.06_wp,                  & !  1
     616                                                 0.19_wp, 0.28_wp, 0.09_wp,                  & !  2
     617                                                 0.23_wp, 0.33_wp, 0.11_wp,                  & !  3
     618                                                 0.23_wp, 0.33_wp, 0.11_wp,                  & !  4
     619                                                 0.25_wp, 0.34_wp, 0.14_wp,                  & !  5
     620                                                 0.14_wp, 0.22_wp, 0.06_wp,                  & !  6
     621                                                 0.17_wp, 0.27_wp, 0.06_wp,                  & !  7
     622                                                 0.19_wp, 0.31_wp, 0.06_wp,                  & !  8
     623                                                 0.14_wp, 0.22_wp, 0.06_wp,                  & !  9
     624                                                 0.18_wp, 0.28_wp, 0.06_wp,                  & ! 10
     625                                                 0.43_wp, 0.51_wp, 0.35_wp,                  & ! 11
     626                                                 0.32_wp, 0.40_wp, 0.24_wp,                  & ! 12
     627                                                 0.19_wp, 0.27_wp, 0.10_wp,                  & ! 13
     628                                                 0.77_wp, 0.65_wp, 0.90_wp,                  & ! 14
     629                                                 0.77_wp, 0.65_wp, 0.90_wp,                  & ! 15
     630                                                 0.82_wp, 0.70_wp, 0.95_wp,                  & ! 16
     631                                                 0.08_wp, 0.08_wp, 0.08_wp,                  & ! 17
     632                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 18
     633                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 19
     634                                                 0.30_wp, 0.30_wp, 0.30_wp,                  & ! 20
     635                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 21
     636                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 22
     637                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 23
     638                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 24
     639                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 25
     640                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 26
     641                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 27
     642                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 28
     643                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 29
     644                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 30
     645                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 31
     646                                                 0.17_wp, 0.17_wp, 0.17_wp,                  & ! 32
     647                                                 0.17_wp, 0.17_wp, 0.17_wp                   & ! 33
     648                                               /), (/ 3, 33 /) )
     649
     650    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  &
     651                        rad_lw_cs_hr,                   & !< longwave clear sky radiation heating rate (K/s)
     652                        rad_lw_cs_hr_av,                & !< average of rad_lw_cs_hr
     653                        rad_lw_hr,                      & !< longwave radiation heating rate (K/s)
     654                        rad_lw_hr_av,                   & !< average of rad_sw_hr
     655                        rad_lw_in,                      & !< incoming longwave radiation (W/m2)
     656                        rad_lw_in_av,                   & !< average of rad_lw_in
     657                        rad_lw_out,                     & !< outgoing longwave radiation (W/m2)
     658                        rad_lw_out_av,                  & !< average of rad_lw_out
     659                        rad_sw_cs_hr,                   & !< shortwave clear sky radiation heating rate (K/s)
     660                        rad_sw_cs_hr_av,                & !< average of rad_sw_cs_hr
     661                        rad_sw_hr,                      & !< shortwave radiation heating rate (K/s)
     662                        rad_sw_hr_av,                   & !< average of rad_sw_hr
     663                        rad_sw_in,                      & !< incoming shortwave radiation (W/m2)
     664                        rad_sw_in_av,                   & !< average of rad_sw_in
     665                        rad_sw_out,                     & !< outgoing shortwave radiation (W/m2)
     666                        rad_sw_out_av                      !< average of rad_sw_out
    586667
    587668
     
    589670!-- Variables and parameters used in RRTMG only
    590671#if defined ( __rrtmg )
    591     CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
     672    CHARACTER(LEN=12) ::  rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
    592673
    593674
    594675!
    595676!-- Flag parameters to be passed to RRTMG (should not be changed until ice phase in clouds is allowed)
    596     INTEGER(iwp), PARAMETER :: rrtm_idrv     = 1, & !< flag for longwave upward flux calculation option (0,1)
    597                                rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
    598                                rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
    599                                rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
    600                                rrtm_inflgsw  = 2, & !< flag for sw cloud optical properties (0,1,2)
    601                                rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
    602                                rrtm_liqflgsw = 1    !< flag for sw liquid droplet specifications
    603 
    604 !
    605 !-- The following variables should be only changed with care, as this will
    606 !-- require further setting of some variables, which is currently not
    607 !-- implemented (aerosols, ice phase).
    608     INTEGER(iwp) :: nzt_rad,           & !< upper vertical limit for radiation calculations
    609                     rrtm_icld = 0,     & !< cloud flag (0: clear sky column, 1: cloudy column)
    610                     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)
    611 
    612     INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
    613 
    614     LOGICAL :: snd_exists = .FALSE.      !< flag parameter to check whether a user-defined input files exists
    615     LOGICAL :: sw_exists = .FALSE.       !< flag parameter to check whether that required rrtmg sw file exists
    616     LOGICAL :: lw_exists = .FALSE.       !< flag parameter to check whether that required rrtmg lw file exists
    617 
    618 
    619     REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
    620 
    621     REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !< hypostatic pressure from sounding data (hPa)
    622                                            rrtm_tsfc,   & !< dummy array for storing surface temperature
    623                                            t_snd          !< actual temperature from sounding data (hPa)
    624 
    625     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rrtm_ccl4vmr,   & !< CCL4 volume mixing ratio (g/mol)
    626                                              rrtm_cfc11vmr,  & !< CFC11 volume mixing ratio (g/mol)
    627                                              rrtm_cfc12vmr,  & !< CFC12 volume mixing ratio (g/mol)
    628                                              rrtm_cfc22vmr,  & !< CFC22 volume mixing ratio (g/mol)
    629                                              rrtm_ch4vmr,    & !< CH4 volume mixing ratio
    630                                              rrtm_cicewp,    & !< in-cloud ice water path (g/m2)
    631                                              rrtm_cldfr,     & !< cloud fraction (0,1)
    632                                              rrtm_cliqwp,    & !< in-cloud liquid water path (g/m2)
    633                                              rrtm_co2vmr,    & !< CO2 volume mixing ratio (g/mol)
    634                                              rrtm_emis,      & !< surface emissivity (0-1)
    635                                              rrtm_h2ovmr,    & !< H2O volume mixing ratio
    636                                              rrtm_n2ovmr,    & !< N2O volume mixing ratio
    637                                              rrtm_o2vmr,     & !< O2 volume mixing ratio
    638                                              rrtm_o3vmr,     & !< O3 volume mixing ratio
    639                                              rrtm_play,      & !< pressure layers (hPa, zu-grid)
    640                                              rrtm_plev,      & !< pressure layers (hPa, zw-grid)
    641                                              rrtm_reice,     & !< cloud ice effective radius (microns)
    642                                              rrtm_reliq,     & !< cloud water drop effective radius (microns)
    643                                              rrtm_tlay,      & !< actual temperature (K, zu-grid)
    644                                              rrtm_tlev,      & !< actual temperature (K, zw-grid)
    645                                              rrtm_lwdflx,    & !< RRTM output of incoming longwave radiation flux (W/m2)
    646                                              rrtm_lwdflxc,   & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
    647                                              rrtm_lwuflx,    & !< RRTM output of outgoing longwave radiation flux (W/m2)
    648                                              rrtm_lwuflxc,   & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
    649                                              rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
    650                                              rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
    651                                              rrtm_lwhr,      & !< RRTM output of longwave radiation heating rate (K/d)
    652                                              rrtm_lwhrc,     & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
    653                                              rrtm_swdflx,    & !< RRTM output of incoming shortwave radiation flux (W/m2)
    654                                              rrtm_swdflxc,   & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
    655                                              rrtm_swuflx,    & !< RRTM output of outgoing shortwave radiation flux (W/m2)
    656                                              rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
    657                                              rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
    658                                              rrtm_swhrc,     & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
    659                                              rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m2)
    660                                              rrtm_difdflux     !< RRTM output of incoming diffuse shortwave (W/m2)
    661 
    662     REAL(wp), DIMENSION(1) ::                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
    663                                              rrtm_aldir,     & !< surface albedo for longwave direct radiation
    664                                              rrtm_asdif,     & !< surface albedo for shortwave diffuse radiation
    665                                              rrtm_asdir        !< surface albedo for shortwave direct radiation
     677    INTEGER(iwp), PARAMETER ::  rrtm_idrv     = 1, & !< flag for longwave upward flux calculation option (0,1)
     678                                rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
     679                                rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
     680                                rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
     681                                rrtm_inflgsw  = 2, & !< flag for sw cloud optical properties (0,1,2)
     682                                rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
     683                                rrtm_liqflgsw = 1     !< flag for sw liquid droplet specifications
     684
     685!
     686!-- The following variables should be only changed with care, as this will require further setting
     687!-- of some variables, which is currently not implemented (aerosols, ice phase).
     688    INTEGER(iwp) :: nzt_rad,       &  !< upper vertical limit for radiation calculations
     689                    rrtm_icld = 0, &  !< cloud flag (0: clear sky column, 1: cloudy column)
     690                    rrtm_iaer = 0     !< aerosol option flag (0: no aerosol layers, for lw only: 6
     691                                      !< (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented)
     692    INTEGER(iwp) ::  nc_stat          !< local variable for storin the result of netCDF calls for error message handling
     693
     694    LOGICAL ::  snd_exists = .FALSE.  !< flag parameter to check whether a user-defined input files exists
     695    LOGICAL ::  sw_exists = .FALSE.   !< flag parameter to check whether that required rrtmg sw file exists
     696    LOGICAL ::  lw_exists = .FALSE.   !< flag parameter to check whether that required rrtmg lw file exists
     697
     698
     699    REAL(wp), PARAMETER ::  mol_mass_air_d_wv = 1.607793_wp  !< molecular weight dry air / water vapor
     700
     701    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  hyp_snd,   &       !< hypostatic pressure from sounding data (hPa)
     702                                              rrtm_tsfc, &       !< dummy array for storing surface temperature
     703                                              t_snd              !< actual temperature from sounding data (hPa)
     704    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rrtm_ccl4vmr,   &  !< CCL4 volume mixing ratio (g/mol)
     705                                              rrtm_cfc11vmr,  &  !< CFC11 volume mixing ratio (g/mol)
     706                                              rrtm_cfc12vmr,  &  !< CFC12 volume mixing ratio (g/mol)
     707                                              rrtm_cfc22vmr,  &  !< CFC22 volume mixing ratio (g/mol)
     708                                              rrtm_ch4vmr,    &  !< CH4 volume mixing ratio
     709                                              rrtm_cicewp,    &  !< in-cloud ice water path (g/m2)
     710                                              rrtm_cldfr,     &  !< cloud fraction (0,1)
     711                                              rrtm_cliqwp,    &  !< in-cloud liquid water path (g/m2)
     712                                              rrtm_co2vmr,    &  !< CO2 volume mixing ratio (g/mol)
     713                                              rrtm_emis,      &  !< surface emissivity (0-1)
     714                                              rrtm_h2ovmr,    &  !< H2O volume mixing ratio
     715                                              rrtm_n2ovmr,    &  !< N2O volume mixing ratio
     716                                              rrtm_o2vmr,     &  !< O2 volume mixing ratio
     717                                              rrtm_o3vmr,     &  !< O3 volume mixing ratio
     718                                              rrtm_play,      &  !< pressure layers (hPa, zu-grid)
     719                                              rrtm_plev,      &  !< pressure layers (hPa, zw-grid)
     720                                              rrtm_reice,     &  !< cloud ice effective radius (microns)
     721                                              rrtm_reliq,     &  !< cloud water drop effective radius (microns)
     722                                              rrtm_tlay,      &  !< actual temperature (K, zu-grid)
     723                                              rrtm_tlev,      &  !< actual temperature (K, zw-grid)
     724                                              rrtm_lwdflx,    &  !< RRTM output of incoming longwave radiation flux (W/m2)
     725                                              rrtm_lwdflxc,   &  !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
     726                                              rrtm_lwuflx,    &  !< RRTM output of outgoing longwave radiation flux (W/m2)
     727                                              rrtm_lwuflxc,   &  !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
     728                                              rrtm_lwuflx_dt, &  !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
     729                                              rrtm_lwuflxc_dt,&  !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
     730                                              rrtm_lwhr,      &  !< RRTM output of longwave radiation heating rate (K/d)
     731                                              rrtm_lwhrc,     &  !< RRTM output of incoming longwave clear sky radiation heating
     732                                                                 !< rate (K/d)
     733                                              rrtm_swdflx,    &  !< RRTM output of incoming shortwave radiation flux (W/m2)
     734                                              rrtm_swdflxc,   &  !< RRTM output of outgoing clear sky shortwave radiation flux(W/m2)
     735                                              rrtm_swuflx,    &  !< RRTM output of outgoing shortwave radiation flux (W/m2)
     736                                              rrtm_swuflxc,   &  !< RRTM output of incoming clear sky shortwave radiation flux(W/m2)
     737                                              rrtm_swhr,      &  !< RRTM output of shortwave radiation heating rate (K/d)
     738                                              rrtm_swhrc,     &  !< RRTM output of incoming shortwave clear sky radiation heating
     739                                                                 !< rate (K/d)
     740                                              rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m2)
     741                                              rrtm_difdflux      !< RRTM output of incoming diffuse shortwave (W/m2)
     742
     743    REAL(wp), DIMENSION(1) ::  rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
     744                               rrtm_aldir,     & !< surface albedo for longwave direct radiation
     745                               rrtm_asdif,     & !< surface albedo for shortwave diffuse radiation
     746                               rrtm_asdir         !< surface albedo for shortwave direct radiation
    666747
    667748!
    668749!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
    669     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   & !< incoming clear sky longwave radiation (W/m2) (not used)
    670                                                 rad_lw_cs_out,  & !< outgoing clear sky longwave radiation (W/m2) (not used)
    671                                                 rad_sw_cs_in,   & !< incoming clear sky shortwave radiation (W/m2) (not used)
    672                                                 rad_sw_cs_out,  & !< outgoing clear sky shortwave radiation (W/m2) (not used)
    673                                                 rrtm_lw_tauaer, & !< lw aerosol optical depth
    674                                                 rrtm_lw_taucld, & !< lw in-cloud optical depth
    675                                                 rrtm_sw_taucld, & !< sw in-cloud optical depth
    676                                                 rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
    677                                                 rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
    678                                                 rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
    679                                                 rrtm_sw_tauaer, & !< sw aerosol optical depth
    680                                                 rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
    681                                                 rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
    682                                                 rrtm_sw_ecaer     !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
     750    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   &  !< incoming clear sky longwave radiation (W/m2) (not used)
     751                                                rad_lw_cs_out,  &  !< outgoing clear sky longwave radiation (W/m2) (not used)
     752                                                rad_sw_cs_in,   &  !< incoming clear sky shortwave radiation (W/m2) (not used)
     753                                                rad_sw_cs_out,  &  !< outgoing clear sky shortwave radiation (W/m2) (not used)
     754                                                rrtm_lw_tauaer, &  !< lw aerosol optical depth
     755                                                rrtm_lw_taucld, &  !< lw in-cloud optical depth
     756                                                rrtm_sw_taucld, &  !< sw in-cloud optical depth
     757                                                rrtm_sw_ssacld, &  !< sw in-cloud single scattering albedo
     758                                                rrtm_sw_asmcld, &  !< sw in-cloud asymmetry parameter
     759                                                rrtm_sw_fsfcld, &  !< sw in-cloud forward scattering fraction
     760                                                rrtm_sw_tauaer, &  !< sw aerosol optical depth
     761                                                rrtm_sw_ssaaer, &  !< sw aerosol single scattering albedo
     762                                                rrtm_sw_asmaer, &  !< sw aerosol asymmetry parameter
     763                                                rrtm_sw_ecaer      !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
    683764
    684765#endif
    685766!
    686767!-- Parameters of urban and land surface models
    687     INTEGER(iwp)                                   ::  nz_urban                           !< number of layers of urban surface (will be calculated)
    688     INTEGER(iwp)                                   ::  nz_plant                           !< number of layers of plant canopy (will be calculated)
    689     INTEGER(iwp)                                   ::  nz_urban_b                         !< bottom layer of urban surface (will be calculated)
    690     INTEGER(iwp)                                   ::  nz_urban_t                         !< top layer of urban surface (will be calculated)
    691     INTEGER(iwp)                                   ::  nz_plant_t                         !< top layer of plant canopy (will be calculated)
    692 !-- parameters of urban and land surface models
    693     INTEGER(iwp), PARAMETER                        ::  nzut_free = 3                      !< number of free layers above top of of topography
    694     INTEGER(iwp), PARAMETER                        ::  ndsvf = 2                          !< number of dimensions of real values in SVF
    695     INTEGER(iwp), PARAMETER                        ::  idsvf = 2                          !< number of dimensions of integer values in SVF
    696     INTEGER(iwp), PARAMETER                        ::  ndcsf = 1                          !< number of dimensions of real values in CSF
    697     INTEGER(iwp), PARAMETER                        ::  idcsf = 2                          !< number of dimensions of integer values in CSF
    698     INTEGER(iwp), PARAMETER                        ::  kdcsf = 4                          !< number of dimensions of integer values in CSF calculation array
    699     INTEGER(iwp), PARAMETER                        ::  id = 1                             !< position of d-index in surfl and surf
    700     INTEGER(iwp), PARAMETER                        ::  iz = 2                             !< position of k-index in surfl and surf
    701     INTEGER(iwp), PARAMETER                        ::  iy = 3                             !< position of j-index in surfl and surf
    702     INTEGER(iwp), PARAMETER                        ::  ix = 4                             !< position of i-index in surfl and surf
    703     INTEGER(iwp), PARAMETER                        ::  im = 5                             !< position of surface m-index in surfl and surf
    704     INTEGER(iwp), PARAMETER                        ::  nidx_surf = 5                      !< number of indices in surfl and surf
    705 
    706     INTEGER(iwp), PARAMETER                        ::  nsurf_type = 10                    !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
    707 
    708     INTEGER(iwp), PARAMETER                        ::  iup_u    = 0                       !< 0 - index of urban upward surface (ground or roof)
    709     INTEGER(iwp), PARAMETER                        ::  idown_u  = 1                       !< 1 - index of urban downward surface (overhanging)
    710     INTEGER(iwp), PARAMETER                        ::  inorth_u = 2                       !< 2 - index of urban northward facing wall
    711     INTEGER(iwp), PARAMETER                        ::  isouth_u = 3                       !< 3 - index of urban southward facing wall
    712     INTEGER(iwp), PARAMETER                        ::  ieast_u  = 4                       !< 4 - index of urban eastward facing wall
    713     INTEGER(iwp), PARAMETER                        ::  iwest_u  = 5                       !< 5 - index of urban westward facing wall
    714 
    715     INTEGER(iwp), PARAMETER                        ::  iup_l    = 6                       !< 6 - index of land upward surface (ground or roof)
    716     INTEGER(iwp), PARAMETER                        ::  inorth_l = 7                       !< 7 - index of land northward facing wall
    717     INTEGER(iwp), PARAMETER                        ::  isouth_l = 8                       !< 8 - index of land southward facing wall
    718     INTEGER(iwp), PARAMETER                        ::  ieast_l  = 9                       !< 9 - index of land eastward facing wall
    719     INTEGER(iwp), PARAMETER                        ::  iwest_l  = 10                      !< 10- index of land westward facing wall
    720 
    721     INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1/)   !< surface normal direction x indices
    722     INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/)   !< surface normal direction y indices
    723     INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/)   !< surface normal direction z indices
    724     REAL(wp),     DIMENSION(0:nsurf_type)          :: facearea                            !< area of single face in respective
    725                                                                                           !< direction (will be calc'd)
    726 
    727 
    728 !-- indices and sizes of urban and land surface models
    729     INTEGER(iwp)                                   ::  startland        !< start index of block of land and roof surfaces
    730     INTEGER(iwp)                                   ::  endland          !< end index of block of land and roof surfaces
    731     INTEGER(iwp)                                   ::  nlands           !< number of land and roof surfaces in local processor
    732     INTEGER(iwp)                                   ::  startwall        !< start index of block of wall surfaces
    733     INTEGER(iwp)                                   ::  endwall          !< end index of block of wall surfaces
    734     INTEGER(iwp)                                   ::  nwalls           !< number of wall surfaces in local processor
    735 
    736 !-- indices needed for RTM netcdf output subroutines
    737     INTEGER(iwp), PARAMETER                        :: nd = 5
    738     CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
    739     INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_u = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
    740     INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_l = (/ iup_l, isouth_l, inorth_l, iwest_l, ieast_l /)
    741     INTEGER(iwp), DIMENSION(0:nd-1)                :: dirstart
    742     INTEGER(iwp), DIMENSION(0:nd-1)                :: dirend
    743 
    744 !-- indices and sizes of urban and land surface models
    745     INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surfl            !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x, m]
    746     INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_linear     !< dtto (linearly allocated array)
    747     INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surf             !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x, m]
    748     INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surf_linear      !< dtto (linearly allocated array)
    749     INTEGER(iwp)                                   ::  nsurfl           !< number of all surfaces in local processor
    750     INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  nsurfs           !< array of number of all surfaces in individual processors
    751     INTEGER(iwp)                                   ::  nsurf            !< global number of surfaces in index array of surfaces (nsurf = proc nsurfs)
    752     INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfstart        !< starts of blocks of surfaces for individual processors in array surf (indexed from 1)
    753                                                                         !< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1]
    754 
    755 !-- block variables needed for calculation of the plant canopy model inside the urban surface model
    756     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pct              !< top layer of the plant canopy
    757     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pch              !< heights of the plant canopy
    758     INTEGER(iwp)                                   ::  npcbl = 0        !< number of the plant canopy gridboxes in local processor
    759     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pcbl             !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i]
    760     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw          !< array of absorbed sw radiation for local plant canopy box
    761     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir       !< array of absorbed direct sw radiation for local plant canopy box
    762     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif       !< array of absorbed diffusion sw radiation for local plant canopy box
    763     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw          !< array of absorbed lw radiation for local plant canopy box
    764 
    765 !-- configuration parameters (they can be setup in PALM config)
    766     LOGICAL                                        ::  raytrace_mpi_rma = .TRUE.          !< use MPI RMA to access LAD and gridsurf from remote processes during raytracing
    767     LOGICAL                                        ::  rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for
    768                                                                                           !< reflected radiation (as opposed to all mutually visible pairs)
    769     LOGICAL                                        ::  plant_lw_interact = .TRUE.         !< whether plant canopy interacts with LW radiation (in addition to SW)
    770     INTEGER(iwp)                                   ::  mrt_nlevels = 0                    !< number of vertical boxes above surface for which to calculate MRT
    771     LOGICAL                                        ::  mrt_skip_roof = .TRUE.             !< do not calculate MRT above roof surfaces
    772     LOGICAL                                        ::  mrt_include_sw = .TRUE.            !< should MRT calculation include SW radiation as well?
    773     INTEGER(wp)                                    ::  mrt_geom = 1                       !< method for MRT direction weights simulating a sphere or a human body
    774     REAL(wp), DIMENSION(2)                         ::  mrt_geom_params = (/ .12_wp, .88_wp /)   !< parameters for the selected method
    775     INTEGER(iwp)                                   ::  nrefsteps = 3                      !< number of reflection steps to perform
    776     REAL(wp), PARAMETER                            ::  ext_coef = 0.6_wp                  !< extinction coefficient (a.k.a. alpha)
    777     INTEGER(iwp), PARAMETER                        ::  rad_version_len = 10               !< length of identification string of rad version
    778     CHARACTER(rad_version_len), PARAMETER          ::  rad_version = 'RAD v. 3.0'         !< identification of version of binary svf and restart files
    779     INTEGER(iwp)                                   ::  raytrace_discrete_elevs = 40       !< number of discretization steps for elevation (nadir to zenith)
    780     INTEGER(iwp)                                   ::  raytrace_discrete_azims = 80       !< number of discretization steps for azimuth (out of 360 degrees)
    781     REAL(wp)                                       ::  max_raytracing_dist = -999.0_wp    !< maximum distance for raytracing (in metres)
    782     REAL(wp)                                       ::  min_irrf_value = 1e-6_wp           !< minimum potential irradiance factor value for raytracing
    783     REAL(wp), DIMENSION(1:30)                      ::  svfnorm_report_thresh = 1e21_wp    !< thresholds of SVF normalization values to report
    784     INTEGER(iwp)                                   ::  svfnorm_report_num                 !< number of SVF normalization thresholds to report
    785 
    786 !-- radiation related arrays to be used in radiation_interaction routine
    787     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rad_sw_in_dir    !< direct sw radiation
    788     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rad_sw_in_diff   !< diffusion sw radiation
    789     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rad_lw_in_diff   !< diffusion lw radiation
    790 
    791 !-- parameters required for RRTMG lower boundary condition
    792     REAL(wp)                   :: albedo_urb      !< albedo value retuned to RRTMG boundary cond.
    793     REAL(wp)                   :: emissivity_urb  !< emissivity value retuned to RRTMG boundary cond.
    794     REAL(wp)                   :: t_rad_urb       !< temperature value retuned to RRTMG boundary cond.
    795 
    796 !-- type for calculation of svf
     768    INTEGER(iwp) ::  nz_urban    !< number of layers of urban surface (will be calculated)
     769    INTEGER(iwp) ::  nz_plant    !< number of layers of plant canopy (will be calculated)
     770    INTEGER(iwp) ::  nz_urban_b  !< bottom layer of urban surface (will be calculated)
     771    INTEGER(iwp) ::  nz_urban_t  !< top layer of urban surface (will be calculated)
     772    INTEGER(iwp) ::  nz_plant_t  !< top layer of plant canopy (will be calculated)
     773!-- Parameters of urban and land surface models
     774    INTEGER(iwp), PARAMETER ::  nzut_free = 3    !< number of free layers above top of of topography
     775    INTEGER(iwp), PARAMETER ::  ndsvf = 2        !< number of dimensions of real values in SVF
     776    INTEGER(iwp), PARAMETER ::  idsvf = 2        !< number of dimensions of integer values in SVF
     777    INTEGER(iwp), PARAMETER ::  ndcsf = 1        !< number of dimensions of real values in CSF
     778    INTEGER(iwp), PARAMETER ::  idcsf = 2        !< number of dimensions of integer values in CSF
     779    INTEGER(iwp), PARAMETER ::  kdcsf = 4        !< number of dimensions of integer values in CSF calculation array
     780    INTEGER(iwp), PARAMETER ::  id = 1           !< position of d-index in surfl and surf
     781    INTEGER(iwp), PARAMETER ::  iz = 2           !< position of k-index in surfl and surf
     782    INTEGER(iwp), PARAMETER ::  iy = 3           !< position of j-index in surfl and surf
     783    INTEGER(iwp), PARAMETER ::  ix = 4           !< position of i-index in surfl and surf
     784    INTEGER(iwp), PARAMETER ::  im = 5           !< position of surface m-index in surfl and surf
     785    INTEGER(iwp), PARAMETER ::  nidx_surf = 5    !< number of indices in surfl and surf
     786    INTEGER(iwp), PARAMETER ::  nsurf_type = 10  !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
     787    INTEGER(iwp), PARAMETER ::  iup_u    = 0     !< 0 - index of urban upward surface (ground or roof)
     788    INTEGER(iwp), PARAMETER ::  idown_u  = 1     !< 1 - index of urban downward surface (overhanging)
     789    INTEGER(iwp), PARAMETER ::  inorth_u = 2     !< 2 - index of urban northward facing wall
     790    INTEGER(iwp), PARAMETER ::  isouth_u = 3     !< 3 - index of urban southward facing wall
     791    INTEGER(iwp), PARAMETER ::  ieast_u  = 4     !< 4 - index of urban eastward facing wall
     792    INTEGER(iwp), PARAMETER ::  iwest_u  = 5     !< 5 - index of urban westward facing wall
     793    INTEGER(iwp), PARAMETER ::  iup_l    = 6     !< 6 - index of land upward surface (ground or roof)
     794    INTEGER(iwp), PARAMETER ::  inorth_l = 7     !< 7 - index of land northward facing wall
     795    INTEGER(iwp), PARAMETER ::  isouth_l = 8     !< 8 - index of land southward facing wall
     796    INTEGER(iwp), PARAMETER ::  ieast_l  = 9     !< 9 - index of land eastward facing wall
     797    INTEGER(iwp), PARAMETER ::  iwest_l  = 10    !< 10- index of land westward facing wall
     798
     799    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1/)  !< surface normal direction x indic.
     800    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/)  !< surface normal direction y indic.
     801    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/)  !< surface normal direction z indic.
     802
     803    REAL(wp), DIMENSION(0:nsurf_type) ::  facearea  !< area of single face in respective direction (will be calc'd)
     804
     805!
     806!-- Indices and sizes of urban and land surface models
     807    INTEGER(iwp) ::  startland  !< start index of block of land and roof surfaces
     808    INTEGER(iwp) ::  endland    !< end index of block of land and roof surfaces
     809    INTEGER(iwp) ::  nlands     !< number of land and roof surfaces in local processor
     810    INTEGER(iwp) ::  startwall  !< start index of block of wall surfaces
     811    INTEGER(iwp) ::  endwall    !< end index of block of wall surfaces
     812    INTEGER(iwp) ::  nwalls     !< number of wall surfaces in local processor
     813!
     814!-- Indices needed for RTM netcdf output subroutines
     815    INTEGER(iwp), PARAMETER                        ::  nd = 5  !<
     816
     817    CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER ::  dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)  !<
     818
     819    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     ::  dirint_u = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)  !<
     820    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     ::  dirint_l = (/ iup_l, isouth_l, inorth_l, iwest_l, ieast_l /)  !<
     821
     822    INTEGER(iwp), DIMENSION(0:nd-1)                ::  dirstart  !<
     823    INTEGER(iwp), DIMENSION(0:nd-1)                ::  dirend    !<
     824!
     825!-- Indices and sizes of urban and land surface models
     826    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surfl         !< coordinates of i-th local surface in local grid
     827                                                                     !< - surfl[:,k] = [d, z, y, x, m]
     828    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_linear  !< dtto (linearly allocated array)
     829    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surf          !< coordinates of i-th surface in grid
     830                                                                     !< - surf[:,k] = [d, z, y, x, m]
     831    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surf_linear   !< dtto (linearly allocated array)
     832    INTEGER(iwp)                                   ::  nsurfl        !< number of all surfaces in local processor
     833    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  nsurfs        !< array of number of all surfaces in individual processors
     834    INTEGER(iwp)                                   ::  nsurf         !< global number of surfaces in index array of surfaces
     835                                                                     !< (nsurf = proc nsurfs)
     836    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfstart     !< starts of blocks of surfaces for individual processors in
     837                                                                     !< array surf (indexed from 1)
     838                          !< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1]
     839!
     840!-- Block variables needed for calculation of the plant canopy model inside the urban surface model
     841    INTEGER(iwp)                              ::  npcbl = 0   !< number of the plant canopy gridboxes in local processor
     842    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pct         !< top layer of the plant canopy
     843    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pch         !< heights of the plant canopy
     844    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pcbl        !< k,j,i coordinates of l-th local plant canopy box
     845                                                              !< pcbl[:,l] = [k, j, i]
     846    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  pcbinsw     !< array of absorbed sw radiation for local plant canopy box
     847    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  pcbinswdir  !< array of absorbed direct sw radiation for local plant canopy box
     848    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  pcbinswdif  !< array of absorbed diffusion sw radiation for local plant canopy box
     849    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  pcbinlw     !< array of absorbed lw radiation for local plant canopy box
     850!
     851!-- Configuration parameters (they can be setup in PALM config)
     852    INTEGER(iwp), PARAMETER ::  rad_version_len = 10  !< length of identification string of rad version
     853
     854    CHARACTER(rad_version_len), PARAMETER ::  rad_version = 'RAD v. 3.0'  !< identification of version of binary svf and restart
     855                                                                          !< files
     856    INTEGER(iwp) ::  mrt_nlevels = 0               !< number of vertical boxes above surface for which to calculate MRT
     857    INTEGER(iwp) ::  svfnorm_report_num            !< number of SVF normalization thresholds to report
     858    INTEGER(iwp) ::  raytrace_discrete_elevs = 40  !< number of discretization steps for elevation (nadir to zenith)
     859    INTEGER(iwp) ::  raytrace_discrete_azims = 80  !< number of discretization steps for azimuth (out of 360 degrees)
     860    INTEGER(wp)  ::  mrt_geom = 1                  !< method for MRT direction weights simulating a sphere or a human body
     861    INTEGER(iwp) ::  nrefsteps = 3                 !< number of reflection steps to perform
     862
     863    LOGICAL ::  raytrace_mpi_rma = .TRUE.            !< use MPI RMA to access LAD and gridsurf from remote processes during
     864                                                     !< raytracing
     865    LOGICAL ::  rad_angular_discretization = .TRUE.  !< whether to use fixed resolution discretization of view factors for
     866                                                     !< reflected radiation (as opposed to all mutually visible pairs)
     867    LOGICAL ::  plant_lw_interact = .TRUE.           !< whether plant canopy interacts with LW radiation (in addition to SW)
     868    LOGICAL ::  mrt_skip_roof = .TRUE.               !< do not calculate MRT above roof surfaces
     869    LOGICAL ::  mrt_include_sw = .TRUE.              !< should MRT calculation include SW radiation as well?
     870
     871    REAL(wp) ::  max_raytracing_dist = -999.0_wp     !< maximum distance for raytracing (in metres)
     872    REAL(wp) ::  min_irrf_value = 1.0E-6_wp          !< minimum potential irradiance factor value for raytracing
     873
     874    REAL(wp), PARAMETER ::  ext_coef = 0.6_wp  !< extinction coefficient (a.k.a. alpha)
     875
     876    REAL(wp), DIMENSION(2)    ::  mrt_geom_params = (/ .12_wp, .88_wp /)  !< parameters for the selected method
     877    REAL(wp), DIMENSION(1:30) ::  svfnorm_report_thresh = 1e21_wp         !< thresholds of SVF normalization values to report
     878!
     879!-- Radiation related arrays to be used in radiation_interaction routine
     880    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_in_dir   !< direct sw radiation
     881    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_in_diff  !< diffusion sw radiation
     882    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_in_diff  !< diffusion lw radiation
     883!
     884!-- Parameters required for RRTMG lower boundary condition
     885    REAL(wp) :: albedo_urb      !< albedo value retuned to RRTMG boundary cond.
     886    REAL(wp) :: emissivity_urb  !< emissivity value retuned to RRTMG boundary cond.
     887    REAL(wp) :: t_rad_urb       !< temperature value retuned to RRTMG boundary cond.
     888!
     889!-- Type for calculation of svf
    797890    TYPE t_svf
    798         INTEGER(iwp)                               :: isurflt           !<
    799         INTEGER(iwp)                               :: isurfs            !<
    800         REAL(wp)                                   :: rsvf              !<
    801         REAL(wp)                                   :: rtransp           !<
     891        INTEGER(iwp) ::  isurflt  !<
     892        INTEGER(iwp) ::  isurfs   !<
     893        REAL(wp)     ::  rsvf     !<
     894        REAL(wp)     ::  rtransp  !<
    802895    END TYPE
    803 
    804 !-- type for calculation of csf
     896!
     897!-- Type for calculation of csf
    805898    TYPE t_csf
    806         INTEGER(iwp)                               :: ip                !<
    807         INTEGER(iwp)                               :: itx               !<
    808         INTEGER(iwp)                               :: ity               !<
    809         INTEGER(iwp)                               :: itz               !<
    810         INTEGER(iwp)                               :: isurfs            !< Idx of source face / -1 for sky
    811         REAL(wp)                                   :: rcvf              !< Canopy view factor for faces /
    812                                                                         !< canopy sink factor for sky (-1)
     899        INTEGER(iwp) ::  ip      !<
     900        INTEGER(iwp) ::  itx     !<
     901        INTEGER(iwp) ::  ity     !<
     902        INTEGER(iwp) ::  itz     !<
     903        INTEGER(iwp) ::  isurfs  !< Idx of source face / -1 for sky
     904        REAL(wp)     ::  rcvf    !< Canopy view factor for faces / canopy sink factor for sky (-1)
    813905    END TYPE
    814 
    815 !-- arrays storing the values of USM
    816     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  svfsurf          !< svfsurf[:,isvf] = index of target and source surface for svf[isvf]
    817     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  svf              !< array of shape view factors+direct irradiation factors for local surfaces
    818     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins          !< array of sw radiation falling to local surface after i-th reflection
    819     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl          !< array of lw radiation for local surface after i-th reflection
    820 
    821     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  skyvf            !< array of sky view factor for each local surface
    822     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  skyvft           !< array of sky view factor including transparency for each local surface
    823     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsitrans         !< dsidir[isvfl,i] = path transmittance of i-th
    824                                                                         !< direction of direct solar irradiance per target surface
    825     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsitransc        !< dtto per plant canopy box
    826     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsidir           !< dsidir[:,i] = unit vector of i-th
    827                                                                         !< direction of direct solar irradiance
    828     INTEGER(iwp)                                   ::  ndsidir          !< number of apparent solar directions used
    829     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  dsidir_rev       !< dsidir_rev[ielev,iazim] = i for dsidir or -1 if not present
    830 
    831     INTEGER(iwp)                                   ::  nmrtbl           !< No. of local grid boxes for which MRT is calculated
    832     INTEGER(iwp)                                   ::  nmrtf            !< number of MRT factors for local processor
    833     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  mrtbl            !< coordinates of i-th local MRT box - surfl[:,i] = [z, y, x]
    834     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  mrtfsurf         !< mrtfsurf[:,imrtf] = index of target MRT box and source surface for mrtf[imrtf]
    835     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtf             !< array of MRT factors for each local MRT box
    836     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtft            !< array of MRT factors including transparency for each local MRT box
    837     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtsky           !< array of sky view factor for each local MRT box
    838     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtskyt          !< array of sky view factor including transparency for each local MRT box
    839     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  mrtdsit          !< array of direct solar transparencies for each local MRT box
    840     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinsw          !< mean SW radiant flux for each MRT box
    841     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinlw          !< mean LW radiant flux for each MRT box
    842     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrt              !< mean radiant temperature for each MRT box
    843     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinsw_av       !< time average mean SW radiant flux for each MRT box
    844     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinlw_av       !< time average mean LW radiant flux for each MRT box
    845     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrt_av           !< time average mean radiant temperature for each MRT box
    846 
    847     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw         !< array of sw radiation falling to local surface including radiation from reflections
    848     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw         !< array of lw radiation falling to local surface including radiation from reflections
    849     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir      !< array of direct sw radiation falling to local surface
    850     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif      !< array of diffuse sw radiation from sky and model boundary falling to local surface
    851     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif      !< array of diffuse lw radiation from sky and model boundary falling to local surface
    852 
    853                                                                         !< Outward radiation is only valid for nonvirtual surfaces
    854     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsl        !< array of reflected sw radiation for local surface in i-th reflection
    855     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutll        !< array of reflected + emitted lw radiation for local surface in i-th reflection
    856     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfouts         !< array of reflected sw radiation for all surfaces in i-th reflection
    857     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutl         !< array of reflected + emitted lw radiation for all surfaces in i-th reflection
    858     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlg         !< global array of incoming lw radiation from plant canopy
    859     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw        !< array of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
    860     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw        !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
    861     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfemitlwl      !< array of emitted lw radiation for local surface used to calculate effective surface temperature for radiation model
    862 
    863 !-- block variables needed for calculation of the plant canopy model inside the urban surface model
    864     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  csfsurf          !< csfsurf[:,icsf] = index of target surface and csf grid index for csf[icsf]
    865     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  csf              !< array of plant canopy sink fators + direct irradiation factors (transparency)
    866     REAL(wp), DIMENSION(:,:,:), POINTER            ::  sub_lad          !< subset of lad_s within urban surface, transformed to plain Z coordinate
     906!
     907!-- Arrays storing the values of USM
     908    INTEGER(iwp)                              ::  ndsidir      !< number of apparent solar directions used
     909    INTEGER(iwp)                              ::  nmrtbl       !< No. of local grid boxes for which MRT is calculated
     910    INTEGER(iwp)                              ::  nmrtf        !< number of MRT factors for local processor
     911    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  mrtbl        !< coordinates of i-th local MRT box - surfl[:,i] = [z, y, x]
     912    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  mrtfsurf     !< mrtfsurf[:,imrtf] = index of target MRT box and source surface for
     913    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  dsidir_rev   !< dsidir_rev[ielev,iazim] = i for dsidir or -1 if not present
     914    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  svfsurf      !< svfsurf[:,isvf] = index of target and source surface for svf[isvf]
     915                                                               !< mrtf[imrtf]
     916
     917
     918    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfins      !< array of sw radiation falling to local surface after i-th
     919                                                               !< reflection
     920    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinl      !< array of lw radiation for local surface after i-th reflection
     921    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  skyvf        !< array of sky view factor for each local surface
     922    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  skyvft       !< array of sky view factor including transparency for each local
     923                                                               !< surface
     924    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtf         !< array of MRT factors for each local MRT box
     925    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtft        !< array of MRT factors including transparency for each local MRT box
     926    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtsky       !< array of sky view factor for each local MRT box
     927    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtskyt      !< array of sky view factor including transparency for each local
     928                                                               !< MRT box
     929    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtinsw      !< mean SW radiant flux for each MRT box
     930    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtinlw      !< mean LW radiant flux for each MRT box
     931    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrt          !< mean radiant temperature for each MRT box
     932    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtinsw_av   !< time average mean SW radiant flux for each MRT box
     933    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrtinlw_av   !< time average mean LW radiant flux for each MRT box
     934    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mrt_av       !< time average mean radiant temperature for each MRT box
     935    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinsw     !< array of sw radiation falling to local surface including radiation
     936                                                               !< from reflections
     937    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinlw     !< array of lw radiation falling to local surface including radiation
     938                                                               !< from reflections
     939    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinswdir  !< array of direct sw radiation falling to local surface
     940    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinswdif  !< array of diffuse sw radiation from sky and model boundary falling
     941                                                               !< to local surface
     942    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinlwdif  !< array of diffuse lw radiation from sky and model boundary falling
     943                                                               !< to local surface
     944                                                               !< Outward radiation is only valid for nonvirtual surfaces
     945    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfoutsl    !< array of reflected sw radiation for local surface in i-th
     946                                                               !< reflection
     947    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfoutll    !< array of reflected + emitted lw radiation for local surface in
     948                                                               !< i-th reflection
     949    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfouts     !< array of reflected sw radiation for all surfaces in i-th
     950                                                               !< reflection
     951    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfoutl     !< array of reflected + emitted lw radiation for all surfaces in
     952                                                               !< i-th reflection
     953    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfinlg     !< global array of incoming lw radiation from plant canopy
     954    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfoutsw    !< array of total sw radiation outgoing from nonvirtual surfaces
     955                                                               !< surfaces after all reflection
     956    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfoutlw    !< array of total lw radiation outgoing from nonvirtual surfaces
     957                                                               !< surfaces after all reflection
     958    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  surfemitlwl  !< array of emitted lw radiation for local surface used to calculate
     959                                                               !< effective surface temperature for radiation model
     960
     961
     962    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  svf          !< array of shape view factors+direct irradiation factors for local
     963                                                               !< surfaces
     964    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  mrtdsit      !< array of direct solar transparencies for each local MRT box
     965    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  dsitrans     !< dsidir[isvfl,i] = path transmittance of i-th
     966                                                               !< direction of direct solar irradiance per target surface
     967    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  dsitransc    !< dtto per plant canopy box
     968    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  dsidir       !< dsidir[:,i] = unit vector of i-th
     969                                                               !< direction of direct solar irradiance
     970!
     971!-- Block variables needed for calculation of the plant canopy model inside the urban surface model
     972    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  csfsurf  !< csfsurf[:,icsf] = index of target surface and csf grid index
     973                                                           !< for csf[icsf]
     974    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  csf      !< array of plant canopy sink fators + direct irradiation factors
     975                                                           !< (transparency)
     976    REAL(wp), DIMENSION(:,:,:), POINTER       ::  sub_lad  !< subset of lad_s within urban surface, transformed to plain
     977                                                           !< Z coordinate
    867978#if defined( __parallel )
    868     REAL(wp), DIMENSION(:), POINTER                ::  sub_lad_g        !< sub_lad globalized (used to avoid MPI RMA calls in raytracing)
     979    REAL(wp), DIMENSION(:), POINTER ::  sub_lad_g  !< sub_lad globalized (used to avoid MPI RMA calls in raytracing)
    869980#endif
    870     REAL(wp)                                       ::  prototype_lad    !< prototype leaf area density for computing effective optical depth
    871     INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nzterr, plantt   !< temporary global arrays for raytracing
    872     INTEGER(iwp)                                   ::  plantt_max
    873 
    874 !-- arrays and variables for calculation of svf and csf
    875     TYPE(t_svf), DIMENSION(:), POINTER             ::  asvf             !< pointer to growing svc array
    876     TYPE(t_csf), DIMENSION(:), POINTER             ::  acsf             !< pointer to growing csf array
    877     TYPE(t_svf), DIMENSION(:), POINTER             ::  amrtf            !< pointer to growing mrtf array
    878     TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  asvf1, asvf2     !< realizations of svf array
    879     TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET ::  acsf1, acsf2     !< realizations of csf array
    880     TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  amrtf1, amrtf2   !< realizations of mftf array
    881     INTEGER(iwp)                                   ::  nsvfla           !< dimmension of array allocated for storage of svf in local processor
    882     INTEGER(iwp)                                   ::  ncsfla           !< dimmension of array allocated for storage of csf in local processor
    883     INTEGER(iwp)                                   ::  nmrtfa           !< dimmension of array allocated for storage of mrt
    884     INTEGER(iwp)                                   ::  msvf, mcsf, mmrtf!< mod for swapping the growing array
    885     INTEGER(iwp), PARAMETER                        ::  gasize = 100000_iwp  !< initial size of growing arrays
    886     REAL(wp), PARAMETER                            ::  grow_factor = 1.4_wp !< growth factor of growing arrays
    887     INTEGER(iwp)                                   ::  nsvfl            !< number of svf for local processor
    888     INTEGER(iwp)                                   ::  ncsfl            !< no. of csf in local processor
    889                                                                         !< needed only during calc_svf but must be here because it is
    890                                                                         !< shared between subroutines calc_svf and raytrace
    891     INTEGER(iwp), DIMENSION(:,:,:,:), POINTER      ::  gridsurf         !< reverse index of local surfl[d,k,j,i] (for case rad_angular_discretization)
    892     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE    ::  gridpcbl         !< reverse index of local pcbl[k,j,i]
    893     INTEGER(iwp), PARAMETER                        ::  nsurf_type_u = 6 !< number of urban surf types (used in gridsurf)
    894 
    895 !-- temporary arrays for calculation of csf in raytracing
    896     INTEGER(iwp)                                   ::  maxboxesg        !< max number of boxes ray can cross in the domain
    897     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  boxes            !< coordinates of gridboxes being crossed by ray
    898     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  crlens           !< array of crossing lengths of ray for particular grid boxes
    899     INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  lad_ip           !< array of numbers of process where lad is stored
     981    INTEGER(iwp)                            ::  plantt_max
     982    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nzterr, plantt  !< temporary global arrays for raytracing
     983
     984    REAL(wp)                                ::  prototype_lad   !< prototype leaf area density for computing effective optical depth
     985!
     986!-- Arrays and variables for calculation of svf and csf
     987    TYPE(t_svf), DIMENSION(:), POINTER             ::  asvf                  !< pointer to growing svc array
     988    TYPE(t_csf), DIMENSION(:), POINTER             ::  acsf                  !< pointer to growing csf array
     989    TYPE(t_svf), DIMENSION(:), POINTER             ::  amrtf                 !< pointer to growing mrtf array
     990    TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  asvf1, asvf2          !< realizations of svf array
     991    TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET ::  acsf1, acsf2          !< realizations of csf array
     992    TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  amrtf1, amrtf2        !< realizations of mftf array
     993
     994    INTEGER(iwp) ::  nsvfla             !< dimmension of array allocated for storage of svf in local processor
     995    INTEGER(iwp) ::  ncsfla             !< dimmension of array allocated for storage of csf in local processor
     996    INTEGER(iwp) ::  nmrtfa             !< dimmension of array allocated for storage of mrt
     997    INTEGER(iwp) ::  msvf, mcsf, mmrtf  !< mod for swapping the growing array
     998    INTEGER(iwp) ::  nsvfl              !< number of svf for local processor
     999    INTEGER(iwp) ::  ncsfl              !< no. of csf in local processor needed only during calc_svf but must be here because
     1000                                        !< it is shared between subroutines calc_svf and raytrace
     1001
     1002    INTEGER(iwp), PARAMETER ::  gasize = 100000_iwp   !< initial size of growing arrays
     1003    INTEGER(iwp), PARAMETER ::  nsurf_type_u = 6      !< number of urban surf types (used in gridsurf)
     1004
     1005    INTEGER(iwp), DIMENSION(:,:,:,:), POINTER   ::  gridsurf              !< reverse index of local surfl[d,k,j,i] (for case
     1006                                                                          !< rad_angular_discretization)
     1007    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  gridpcbl              !< reverse index of local pcbl[k,j,i]
     1008
     1009    REAL(wp), PARAMETER ::  grow_factor = 1.4_wp  !< growth factor of growing arrays
     1010
     1011!
     1012!-- Temporary arrays for calculation of csf in raytracing
     1013    INTEGER(iwp) ::  maxboxesg  !< max number of boxes ray can cross in the domain
     1014
     1015    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  boxes      !< coordinates of gridboxes being crossed by ray
     1016    INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  lad_ip     !< array of numbers of process where lad is stored
     1017
     1018    REAL(wp), DIMENSION(:), ALLOCATABLE ::  crlens     !< array of crossing lengths of ray for particular grid boxes
     1019
    9001020#if defined( __parallel )
    901     INTEGER(kind=MPI_ADDRESS_KIND), &
    902                   DIMENSION(:), ALLOCATABLE        ::  lad_disp         !< array of displaycements of lad in local array of proc lad_ip
    903     INTEGER(iwp)                                   ::  win_lad          !< MPI RMA window for leaf area density
    904     INTEGER(iwp)                                   ::  win_gridsurf     !< MPI RMA window for reverse grid surface index
    905     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  lad_s_ray        !< array of received lad_s for appropriate gridboxes crossed by ray
     1021    INTEGER(iwp) ::  win_lad       !< MPI RMA window for leaf area density
     1022    INTEGER(iwp) ::  win_gridsurf  !< MPI RMA window for reverse grid surface index
     1023
     1024    INTEGER(KIND=MPI_ADDRESS_KIND), DIMENSION(:), ALLOCATABLE ::  lad_disp  !< array of displaycements of lad in local array of
     1025                                                                            !< proc lad_ip
     1026    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad_s_ray  !< array of received lad_s for appropriate gridboxes crossed by ray
    9061027#endif
    907     INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  target_surfl
    908     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  rt2_track
    909     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rt2_track_lad
    910     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rt2_track_dist
    911     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rt2_dist
    912 
    913 !-- arrays for time averages
    914     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfradnet_av    !< average of net radiation to local surface including radiation from reflections
    915     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw_av      !< average of sw radiation falling to local surface including radiation from reflections
    916     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw_av      !< average of lw radiation falling to local surface including radiation from reflections
    917     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir_av   !< average of direct sw radiation falling to local surface
    918     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif_av   !< average of diffuse sw radiation from sky and model boundary falling to local surface
    919     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif_av   !< average of diffuse lw radiation from sky and model boundary falling to local surface
    920     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswref_av   !< average of sw radiation falling to surface from reflections
    921     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwref_av   !< average of lw radiation falling to surface from reflections
    922     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw_av     !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
    923     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw_av     !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
    924     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
    925     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
    926     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw_av       !< Average of pcbinlw
    927     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw_av       !< Average of pcbinsw
    928     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir_av    !< Average of pcbinswdir
    929     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif_av    !< Average of pcbinswdif
    930     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswref_av    !< Average of pcbinswref
     1028    INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::  target_surfl    !<
     1029    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  rt2_track       !<
     1030
     1031    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rt2_track_dist  !<
     1032    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rt2_dist        !<
     1033    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rt2_track_lad   !<
     1034!
     1035!-- Arrays for time averages
     1036    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfradnet_av   !< average of net radiation to local surface including radiation from
     1037                                                            !< reflections
     1038    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinsw_av     !< average of sw radiation falling to local surface including radiation
     1039                                                            !< from reflections
     1040    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlw_av     !< average of lw radiation falling to local surface including radiation
     1041                                                            !< from reflections
     1042    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinswdir_av  !< average of direct sw radiation falling to local surface
     1043    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinswdif_av  !< average of diffuse sw radiation from sky and model boundary falling
     1044                                                            !< to local surface
     1045    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlwdif_av  !< average of diffuse lw radiation from sky and model boundary falling
     1046                                                            !< to local surface
     1047    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinswref_av  !< average of sw radiation falling to surface from reflections
     1048    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlwref_av  !< average of lw radiation falling to surface from reflections
     1049    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutsw_av    !< average of total sw radiation outgoing from nonvirtual surfaces
     1050                                                            !< surfaces after all reflection
     1051    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutlw_av    !< average of total lw radiation outgoing from nonvirtual surfaces
     1052                                                            !< surfaces after all reflection
     1053    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfins_av      !< average of array of residua of sw radiation absorbed in surface after
     1054                                                            !< last reflection
     1055    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinl_av      !< average of array of residua of lw radiation absorbed in surface after
     1056                                                            !< last reflection
     1057    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pcbinlw_av      !< Average of pcbinlw
     1058    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pcbinsw_av      !< Average of pcbinsw
     1059    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pcbinswdir_av   !< Average of pcbinswdir
     1060    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pcbinswdif_av   !< Average of pcbinswdif
     1061    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pcbinswref_av   !< Average of pcbinswref
    9311062
    9321063
     
    9341065!-- Energy balance variables
    9351066!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    936 !-- parameters of the land, roof and wall surfaces
    937     REAL(wp), DIMENSION(:), ALLOCATABLE            :: albedo_surf        !< albedo of the surface
    938     REAL(wp), DIMENSION(:), ALLOCATABLE            :: emiss_surf         !< emissivity of the wall surface
    939 !
    940 !-- External radiation. Depending on the given level of detail either a 1D or
    941 !-- a 3D array will be allocated.
    942     TYPE( real_1d_3d ) ::  rad_lw_in_f     !< external incoming longwave radiation, from observation or model
    943     TYPE( real_1d_3d ) ::  rad_sw_in_f     !< external incoming shortwave radiation, from observation or model
    944     TYPE( real_1d_3d ) ::  rad_sw_in_dif_f !< external incoming shortwave radiation, diffuse part, from observation or model
    945     TYPE( real_1d_3d ) ::  time_rad_f      !< time dimension for external radiation, from observation or model
     1067!-- Parameters of the land, roof and wall surfaces
     1068    REAL(wp), DIMENSION(:), ALLOCATABLE :: albedo_surf  !< albedo of the surface
     1069    REAL(wp), DIMENSION(:), ALLOCATABLE :: emiss_surf   !< emissivity of the wall surface
     1070!
     1071!-- External radiation. Depending on the given level of detail either a 1D or a 3D array will be
     1072!-- allocated.
     1073    TYPE( real_1d_3d ) ::  rad_lw_in_f      !< external incoming longwave radiation, from observation or model
     1074    TYPE( real_1d_3d ) ::  rad_sw_in_f      !< external incoming shortwave radiation, from observation or model
     1075    TYPE( real_1d_3d ) ::  rad_sw_in_dif_f  !< external incoming shortwave radiation, diffuse part, from observation or model
     1076    TYPE( real_1d_3d ) ::  time_rad_f       !< time dimension for external radiation, from observation or model
    9461077
    9471078    INTERFACE radiation_check_data_output
     
    10561187!
    10571188!-- Public functions / NEEDS SORTING
    1058     PUBLIC radiation_check_data_output, radiation_check_data_output_pr,        &
    1059            radiation_check_data_output_ts,                                     &
    1060            radiation_check_parameters, radiation_control,                      &
    1061            radiation_header, radiation_init, radiation_parin,                  &
    1062            radiation_3d_data_averaging,                                        &
    1063            radiation_data_output_2d, radiation_data_output_3d,                 &
    1064            radiation_define_netcdf_grid, radiation_wrd_local,                  &
    1065            radiation_rrd_local, radiation_data_output_mask,                    &
    1066            radiation_calc_svf, radiation_write_svf,                            &
    1067            radiation_interaction, radiation_interaction_init,                  &
    1068            radiation_read_svf, radiation_presimulate_solar_pos
     1189    PUBLIC radiation_check_data_output,                                                            &
     1190           radiation_check_data_output_pr,                                                         &
     1191           radiation_check_data_output_ts,                                                         &
     1192           radiation_check_parameters,                                                             &
     1193           radiation_control,                                                                      &
     1194           radiation_header,                                                                       &
     1195           radiation_init,                                                                         &
     1196           radiation_parin,                                                                        &
     1197           radiation_3d_data_averaging,                                                            &
     1198           radiation_data_output_2d,                                                               &
     1199           radiation_data_output_3d,                                                               &
     1200           radiation_define_netcdf_grid,                                                           &
     1201           radiation_wrd_local,                                                                    &
     1202           radiation_rrd_local,                                                                    &
     1203           radiation_data_output_mask,                                                             &
     1204           radiation_calc_svf,                                                                     &
     1205           radiation_write_svf,                                                                    &
     1206           radiation_interaction,                                                                  &
     1207           radiation_interaction_init,                                                             &
     1208           radiation_read_svf,                                                                     &
     1209           radiation_presimulate_solar_pos
    10691210
    10701211
    10711212!
    10721213!-- Public variables and constants / NEEDS SORTING
    1073     PUBLIC albedo, albedo_type, decl_1, decl_2, decl_3, dots_rad, dt_radiation,&
    1074            emissivity, force_radiation_call, lat, lon, mrt_geom,               &
    1075            mrt_geom_params,                                                    &
    1076            mrt_include_sw, mrt_nlevels, mrtbl, mrtinsw, mrtinlw, nmrtbl,       &
    1077            rad_net_av, radiation, radiation_scheme, rad_lw_in,                 &
    1078            rad_lw_in_av, rad_lw_out, rad_lw_out_av,                            &
    1079            rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
    1080            rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
    1081            rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, solar_constant,           &
    1082            skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
    1083            cos_zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,   &
    1084            idir, jdir, kdir, id, iz, iy, ix,                                   &
    1085            iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
    1086            iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
    1087            nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf,                 &
    1088            idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
    1089            radiation_interactions, startwall, startland, endland, endwall,     &
    1090            skyvf, skyvft, radiation_interactions_on, average_radiation,        &
    1091            rad_sw_in_diff, rad_sw_in_dir
     1214    PUBLIC albedo,                                                                                 &
     1215           albedo_type,                                                                            &
     1216           average_radiation,                                                                      &
     1217           calc_zenith,                                                                            &
     1218           cos_zenith,                                                                             &
     1219           decl_1,                                                                                 &
     1220           decl_2,                                                                                 &
     1221           decl_3,                                                                                 &
     1222           dots_rad,                                                                               &
     1223           dt_radiation,                                                                           &
     1224           endland,                                                                                &
     1225           endwall,                                                                                &
     1226           emissivity,                                                                             &
     1227           force_radiation_call,                                                                   &
     1228           id,                                                                                     &
     1229           idcsf,                                                                                  &
     1230           idir,                                                                                   &
     1231           idsvf,                                                                                  &
     1232           ieast_l,                                                                                &
     1233           ieast_u,                                                                                &
     1234           inorth_l,                                                                               &
     1235           inorth_u,                                                                               &
     1236           isouth_l,                                                                               &
     1237           isouth_u,                                                                               &
     1238           iup_l,                                                                                  &
     1239           iup_u,                                                                                  &
     1240           iwest_l,                                                                                &
     1241           iwest_u,                                                                                &
     1242           ix,                                                                                     &
     1243           iy,                                                                                     &
     1244           iz,                                                                                     &
     1245           jdir,                                                                                   &
     1246           kdcsf,                                                                                  &
     1247           kdir,                                                                                   &
     1248           lat,                                                                                    &
     1249           lon,                                                                                    &
     1250           mrtbl,                                                                                  &
     1251           mrt_geom,                                                                               &
     1252           mrt_geom_params,                                                                        &
     1253           mrt_include_sw,                                                                         &
     1254           mrt_nlevels,                                                                            &
     1255           mrtinsw,                                                                                &
     1256           mrtinlw,                                                                                &
     1257           nmrtbl,                                                                                 &
     1258           nsurf_type,                                                                             &
     1259           nz_urban_b,                                                                             &
     1260           nz_urban_t,                                                                             &
     1261           nz_urban,                                                                               &
     1262           nsurf,                                                                                  &
     1263           ndsvf,                                                                                  &
     1264           ndcsf,                                                                                  &
     1265           pch,                                                                                    &
     1266           pct,                                                                                    &
     1267           rad_net_av,                                                                             &
     1268           radiation,                                                                              &
     1269           radiation_scheme,                                                                       &
     1270           rad_lw_in,                                                                              &
     1271           rad_lw_in_av,                                                                           &
     1272           rad_lw_out,                                                                             &
     1273           rad_lw_out_av,                                                                          &
     1274           rad_lw_cs_hr,                                                                           &
     1275           rad_lw_cs_hr_av,                                                                        &
     1276           rad_lw_hr,                                                                              &
     1277           rad_lw_hr_av,                                                                           &
     1278           rad_sw_in,                                                                              &
     1279           rad_sw_in_av,                                                                           &
     1280           rad_sw_out,                                                                             &
     1281           rad_sw_out_av,                                                                          &
     1282           rad_sw_cs_hr,                                                                           &
     1283           rad_sw_cs_hr_av,                                                                        &
     1284           rad_sw_hr,                                                                              &
     1285           rad_sw_hr_av,                                                                           &
     1286           radiation_interactions,                                                                 &
     1287           radiation_interactions_on,                                                              &
     1288           rad_sw_in_diff,                                                                         &
     1289           rad_sw_in_dir,                                                                          &
     1290           solar_constant,                                                                         &
     1291           skip_time_do_radiation,                                                                 &
     1292           sun_direction,                                                                          &
     1293           sun_dir_lat,                                                                            &
     1294           sun_dir_lon,                                                                            &
     1295           startwall,                                                                              &
     1296           startland,                                                                              &
     1297           skyvf,                                                                                  &
     1298           skyvft,                                                                                 &
     1299           time_radiation,                                                                         &
     1300           unscheduled_radiation_calls
     1301
    10921302
    10931303
     
    10991309
    11001310
    1101 !------------------------------------------------------------------------------!
     1311!--------------------------------------------------------------------------------------------------!
    11021312! Description:
    11031313! ------------
    11041314!> This subroutine controls the calls of the radiation schemes
    1105 !------------------------------------------------------------------------------!
    1106     SUBROUTINE radiation_control
    1107 
    1108 
    1109        IMPLICIT NONE
    1110 
    1111 
    1112        IF ( debug_output_timestep )  CALL debug_message( 'radiation_control', 'start' )
    1113 
    1114 
    1115        SELECT CASE ( TRIM( radiation_scheme ) )
    1116 
    1117           CASE ( 'constant' )
    1118              CALL radiation_constant
    1119 
    1120           CASE ( 'clear-sky' )
     1315!--------------------------------------------------------------------------------------------------!
     1316 SUBROUTINE radiation_control
     1317
     1318
     1319    IMPLICIT NONE
     1320
     1321
     1322    IF ( debug_output_timestep )  CALL debug_message( 'radiation_control', 'start' )
     1323
     1324
     1325    SELECT CASE ( TRIM( radiation_scheme ) )
     1326
     1327       CASE ( 'constant' )
     1328          CALL radiation_constant
     1329
     1330       CASE ( 'clear-sky' )
     1331          CALL radiation_clearsky
     1332
     1333       CASE ( 'rrtmg' )
     1334          CALL radiation_rrtmg
     1335
     1336       CASE ( 'external' )
     1337!
     1338!--       During spinup apply clear-sky model
     1339          IF ( time_since_reference_point < 0.0_wp )  THEN
    11211340             CALL radiation_clearsky
    1122 
    1123           CASE ( 'rrtmg' )
    1124              CALL radiation_rrtmg
    1125 
    1126           CASE ( 'external' )
    1127 !
    1128 !--          During spinup apply clear-sky model
    1129              IF ( time_since_reference_point < 0.0_wp )  THEN
    1130                 CALL radiation_clearsky
    1131              ELSE
    1132                 CALL radiation_external
    1133              ENDIF
    1134 
    1135           CASE DEFAULT
    1136 
    1137        END SELECT
    1138 
    1139        IF ( debug_output_timestep )  CALL debug_message( 'radiation_control', 'end' )
    1140 
    1141     END SUBROUTINE radiation_control
    1142 
    1143 !------------------------------------------------------------------------------!
     1341          ELSE
     1342             CALL radiation_external
     1343          ENDIF
     1344
     1345       CASE DEFAULT
     1346
     1347    END SELECT
     1348
     1349    IF ( debug_output_timestep )  CALL debug_message( 'radiation_control', 'end' )
     1350
     1351 END SUBROUTINE radiation_control
     1352
     1353!--------------------------------------------------------------------------------------------------!
    11441354! Description:
    11451355! ------------
    11461356!> Check data output for radiation model
    1147 !------------------------------------------------------------------------------!
    1148     SUBROUTINE radiation_check_data_output( variable, unit, i, ilen, k )
    1149 
    1150 
    1151        USE control_parameters,                                                 &
    1152            ONLY: data_output, message_string
    1153 
    1154        IMPLICIT NONE
    1155 
    1156        LOGICAL                      ::  directional
    1157        CHARACTER(LEN=*)             ::  unit        !<
    1158        CHARACTER(LEN=*)             ::  variable    !<
    1159        CHARACTER(LEN=varnamelength) ::  var         !< TRIM(variable)
    1160        INTEGER(iwp)                 ::  i, k
    1161        INTEGER(iwp)                 ::  ilast_word
    1162        INTEGER(iwp)                 ::  ilen
    1163 
    1164        var = TRIM(variable)
    1165 !
    1166 !--    Identify directional variables
    1167        ilast_word = SCAN(var, '_', back=.TRUE.)
    1168        IF ( ilast_word > 0 )  THEN
    1169           SELECT CASE ( var(ilast_word:) )
    1170              CASE ( '_roof', '_south', '_north', '_west', '_east' )
    1171                 directional = .TRUE.
    1172                 write(9, *) 'vardir', var
    1173                 flush(9)
    1174                 var = var(1:ilast_word-1)
    1175              CASE DEFAULT
    1176                 directional = .FALSE.
    1177                 write(9, *) 'varnd', var
    1178                 flush(9)
    1179           END SELECT
    1180        ELSE
    1181           directional = .FALSE.
    1182        END IF
    1183 
    1184        IF ( directional )  THEN
    1185           SELECT CASE ( var )
    1186              CASE ( 'rtm_rad_net', 'rtm_rad_insw', 'rtm_rad_inlw',             &
    1187                     'rtm_rad_inswdir', 'rtm_rad_inswdif', 'rtm_rad_inswref',   &
    1188                     'rtm_rad_inlwdif', 'rtm_rad_inlwref', 'rtm_rad_outsw',     &
    1189                     'rtm_rad_outlw', 'rtm_rad_ressw', 'rtm_rad_reslw' )
    1190                 IF ( .NOT.  radiation ) THEN
    1191                    message_string = 'output of "' // var // '" require'        &
    1192                                     // 's radiation = .TRUE.'
    1193                    CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
    1194                 ENDIF
    1195                 unit = 'W/m2'
    1196 
    1197              CASE ( 'rtm_svf', 'rtm_dif', 'rtm_skyvf', 'rtm_skyvft',           &
    1198                     'rtm_surfalb', 'rtm_surfemis' )
    1199                 IF ( .NOT.  radiation ) THEN
    1200                    message_string = 'output of "' // var // '" require'&
    1201                                     // 's radiation = .TRUE.'
    1202                    CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
    1203                 ENDIF
    1204                 unit = '1'
    1205 
    1206              CASE DEFAULT
    1207                 unit = 'illegal'
    1208           END SELECT
    1209 
    1210        ELSE
    1211           SELECT CASE ( var )
    1212              CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out',     &
    1213                     'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out'  )
    1214                 IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
    1215                    message_string = '"output of "' // var // '" requi' //       &
    1216                                     'res radiation = .TRUE. and ' //            &
    1217                                     'radiation_scheme = "rrtmg"'
    1218                    CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
    1219                 ENDIF
    1220                 unit = 'K/h'
    1221 
    1222              CASE ( 'rad_net*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',      &
    1223                     'rad_sw_out*' )
    1224                 IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    1225                    message_string = 'illegal value for data_output: "' //       &
    1226                                     var // '" & only 2d-horizontal ' //         &
    1227                                     'cross sections are allowed for this value'
    1228                    CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
    1229                 ENDIF
    1230                 unit = 'W/m2'
    1231 
    1232              CASE ( 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*' )
    1233                 IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    1234                    message_string = 'illegal value for data_output: "' //       &
    1235                                     var // '" & only 2d-horizontal ' //         &
    1236                                     'cross sections are allowed for this value'
    1237                    CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
    1238                 ENDIF
    1239                 IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
    1240                    message_string = 'output of "' // var // '" require'         &
    1241                                     // 's radiation = .TRUE. and radiation_sch' &
    1242                                     // 'eme = "rrtmg"'
    1243                    CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
    1244                 ENDIF
    1245                 unit = ''
    1246 
    1247              CASE ( 'rtm_rad_pc_inlw', 'rtm_rad_pc_insw', 'rtm_rad_pc_inswdir', &
    1248                     'rtm_rad_pc_inswdif', 'rtm_rad_pc_inswref')
    1249                 IF ( .NOT.  radiation ) THEN
    1250                    message_string = 'output of "' // var // '" require'         &
    1251                                     // 's radiation = .TRUE.'
    1252                    CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
    1253                 ENDIF
    1254                 unit = 'W'
    1255 
    1256              CASE ( 'rtm_mrt', 'rtm_mrt_sw', 'rtm_mrt_lw'  )
    1257                 IF ( .NOT.  radiation ) THEN
    1258                    message_string = 'output of "' // var // '" require'         &
    1259                                     // 's radiation = .TRUE.'
    1260                    CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
    1261                 ENDIF
    1262                 IF ( mrt_nlevels == 0 ) THEN
    1263                    message_string = 'output of "' // var // '" require'         &
    1264                                     // 's mrt_nlevels > 0'
    1265                    CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
    1266                 ENDIF
    1267                 IF ( var == 'rtm_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
    1268                    message_string = 'output of "' // var // '" require'         &
    1269                                     // 's rtm_mrt_sw = .TRUE.'
    1270                    CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
    1271                 ENDIF
    1272                 IF ( var == 'rtm_mrt' ) THEN
    1273                    unit = 'K'
    1274                 ELSE
    1275                    unit = 'W m-2'
    1276                 ENDIF
    1277 
    1278              CASE DEFAULT
    1279                 unit = 'illegal'
    1280 
    1281           END SELECT
    1282        END IF
    1283 
    1284     END SUBROUTINE radiation_check_data_output
    1285 
    1286 
    1287 !------------------------------------------------------------------------------!
     1357!--------------------------------------------------------------------------------------------------!
     1358 SUBROUTINE radiation_check_data_output( variable, unit, i, ilen, k )
     1359
     1360
     1361    USE control_parameters,                                                                        &
     1362        ONLY: data_output,                                                                         &
     1363              message_string
     1364
     1365    IMPLICIT NONE
     1366
     1367    CHARACTER(LEN=*)             ::  unit         !<
     1368    CHARACTER(LEN=*)             ::  variable     !<
     1369    CHARACTER(LEN=varnamelength) ::  var          !< TRIM(variable)
     1370
     1371    INTEGER(iwp)                 ::  i, k         !<
     1372    INTEGER(iwp)                 ::  ilast_word   !<
     1373    INTEGER(iwp)                 ::  ilen         !<
     1374
     1375    LOGICAL                      ::  directional  !<
     1376
     1377    var = TRIM( variable )
     1378!
     1379!-- Identify directional variables
     1380    ilast_word = SCAN( var, '_', back = .TRUE. )
     1381    IF ( ilast_word > 0 )  THEN
     1382       SELECT CASE ( var(ilast_word:) )
     1383          CASE ( '_roof', '_south', '_north', '_west', '_east' )
     1384             directional = .TRUE.
     1385             WRITE( 9, * ) 'vardir', var
     1386             FLUSH( 9 )
     1387             var = var(1:ilast_word-1)
     1388          CASE DEFAULT
     1389             directional = .FALSE.
     1390             WRITE( 9, * ) 'varnd', var
     1391             FLUSH( 9 )
     1392       END SELECT
     1393    ELSE
     1394       directional = .FALSE.
     1395    END IF
     1396
     1397    IF ( directional )  THEN
     1398       SELECT CASE ( var )
     1399          CASE ( 'rtm_rad_net', 'rtm_rad_insw', 'rtm_rad_inlw', 'rtm_rad_inswdir',                 &
     1400                 'rtm_rad_inswdif', 'rtm_rad_inswref', 'rtm_rad_inlwdif', 'rtm_rad_inlwref',       &
     1401                 'rtm_rad_outsw', 'rtm_rad_outlw', 'rtm_rad_ressw', 'rtm_rad_reslw' )
     1402             IF ( .NOT.  radiation )  THEN
     1403                message_string = 'output of "' // var // '" requires radiation = .TRUE.'
     1404                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
     1405             ENDIF
     1406             unit = 'W/m2'
     1407
     1408          CASE ( 'rtm_svf', 'rtm_dif', 'rtm_skyvf', 'rtm_skyvft', 'rtm_surfalb', 'rtm_surfemis' )
     1409             IF ( .NOT.  radiation )  THEN
     1410                message_string = 'output of "' // var // '" requires radiation = .TRUE.'
     1411                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
     1412             ENDIF
     1413             unit = '1'
     1414
     1415          CASE DEFAULT
     1416             unit = 'illegal'
     1417       END SELECT
     1418
     1419    ELSE
     1420       SELECT CASE ( var )
     1421          CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out', 'rad_sw_cs_hr',           &
     1422                 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out'  )
     1423             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
     1424                message_string = '"output of "' // var // '" requires radiation = .TRUE. and ' //  &
     1425                                 'radiation_scheme = "rrtmg"'
     1426                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
     1427             ENDIF
     1428             unit = 'K/h'
     1429
     1430          CASE ( 'rad_net*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*', 'rad_sw_out*' )
     1431             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
     1432                message_string = 'illegal value for data_output: "' // var //                      &
     1433                                 '" & only 2d-horizontal cross sections are allowed for this value'
     1434                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
     1435             ENDIF
     1436             unit = 'W/m2'
     1437
     1438          CASE ( 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*' )
     1439             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
     1440                message_string = 'illegal value for data_output: "' // var //                      &
     1441                                 '" & only 2d-horizontal cross sections are allowed for this value'
     1442                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
     1443             ENDIF
     1444             IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
     1445                message_string = 'output of "' // var // '" requires radiation = .TRUE. ' //       &
     1446                                 'and radiation_scheme = "rrtmg"'
     1447                CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
     1448             ENDIF
     1449             unit = ''
     1450
     1451          CASE ( 'rtm_rad_pc_inlw', 'rtm_rad_pc_insw', 'rtm_rad_pc_inswdir', 'rtm_rad_pc_inswdif', &
     1452                 'rtm_rad_pc_inswref')
     1453             IF ( .NOT.  radiation )  THEN
     1454                message_string = 'output of "' // var // '" requires radiation = .TRUE.'
     1455                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
     1456             ENDIF
     1457             unit = 'W'
     1458
     1459          CASE ( 'rtm_mrt', 'rtm_mrt_sw', 'rtm_mrt_lw'  )
     1460             IF ( .NOT.  radiation )  THEN
     1461                message_string = 'output of "' // var // '" requires radiation = .TRUE.'
     1462                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
     1463             ENDIF
     1464             IF ( mrt_nlevels == 0 )  THEN
     1465                message_string = 'output of "' // var // '" requires mrt_nlevels > 0'
     1466                CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
     1467             ENDIF
     1468             IF ( var == 'rtm_mrt_sw'  .AND.  .NOT. mrt_include_sw )  THEN
     1469                message_string = 'output of "' // var // '" requires rtm_mrt_sw = .TRUE.'
     1470                CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
     1471             ENDIF
     1472             IF ( var == 'rtm_mrt' )  THEN
     1473                unit = 'K'
     1474             ELSE
     1475                unit = 'W m-2'
     1476             ENDIF
     1477
     1478          CASE DEFAULT
     1479             unit = 'illegal'
     1480
     1481       END SELECT
     1482    END IF
     1483
     1484 END SUBROUTINE radiation_check_data_output
     1485
     1486
     1487!--------------------------------------------------------------------------------------------------!
    12881488! Description:
    12891489! ------------
    12901490!> Set module-specific timeseries units and labels
    1291 !------------------------------------------------------------------------------!
     1491!--------------------------------------------------------------------------------------------------!
    12921492 SUBROUTINE radiation_check_data_output_ts( dots_max, dots_num )
    12931493
    12941494
    1295     INTEGER(iwp),      INTENT(IN)     ::  dots_max
    1296     INTEGER(iwp),      INTENT(INOUT)  ::  dots_num
     1495    INTEGER(iwp), INTENT(IN)    ::  dots_max  !<
     1496    INTEGER(iwp), INTENT(INOUT) ::  dots_num  !<
    12971497
    12981498!
     
    13011501
    13021502!
    1303 !-- Temporary solution to add LSM and radiation time series to the default
    1304 !-- output
     1503!-- Temporary solution to add LSM and radiation time series to the default output
    13051504    IF ( land_surface  .OR.  radiation )  THEN
    13061505       IF ( TRIM( radiation_scheme ) == 'rrtmg' )  THEN
     
    13141513 END SUBROUTINE radiation_check_data_output_ts
    13151514
    1316 !------------------------------------------------------------------------------!
     1515!--------------------------------------------------------------------------------------------------!
    13171516! Description:
    13181517! ------------
    13191518!> Check data output of profiles for radiation model
    1320 !------------------------------------------------------------------------------!
    1321     SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit,      &
    1322                dopr_unit )
    1323 
    1324        USE arrays_3d,                                                          &
    1325            ONLY: zu
    1326 
    1327        USE control_parameters,                                                 &
    1328            ONLY: data_output_pr, message_string
    1329 
    1330        USE indices
    1331 
    1332        USE profil_parameter
    1333 
    1334        USE statistics
    1335 
    1336        IMPLICIT NONE
    1337 
    1338        CHARACTER (LEN=*) ::  unit      !<
    1339        CHARACTER (LEN=*) ::  variable  !<
    1340        CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
    1341 
    1342        INTEGER(iwp) ::  var_count     !<
    1343 
    1344        SELECT CASE ( TRIM( variable ) )
    1345 
    1346          CASE ( 'rad_net' )
    1347              IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
    1348              THEN
    1349                 message_string = 'data_output_pr = ' //                        &
    1350                                  TRIM( data_output_pr(var_count) ) // ' is' // &
    1351                                  'not available for radiation = .FALSE. or ' //&
    1352                                  'radiation_scheme = "constant"'
    1353                 CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
    1354              ELSE
    1355                 dopr_index(var_count) = 99
    1356                 dopr_unit  = 'W/m2'
    1357                 hom(:,2,99,:)  = SPREAD( zw, 2, statistic_regions+1 )
    1358                 unit = dopr_unit
    1359              ENDIF
    1360 
    1361           CASE ( 'rad_lw_in' )
    1362              IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
    1363              THEN
    1364                 message_string = 'data_output_pr = ' //                        &
    1365                                  TRIM( data_output_pr(var_count) ) // ' is' // &
    1366                                  'not available for radiation = .FALSE. or ' //&
    1367                                  'radiation_scheme = "constant"'
    1368                 CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
    1369              ELSE
    1370                 dopr_index(var_count) = 100
    1371                 dopr_unit  = 'W/m2'
    1372                 hom(:,2,100,:)  = SPREAD( zw, 2, statistic_regions+1 )
    1373                 unit = dopr_unit
    1374              ENDIF
    1375 
    1376           CASE ( 'rad_lw_out' )
    1377              IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
    1378              THEN
    1379                 message_string = 'data_output_pr = ' //                        &
    1380                                  TRIM( data_output_pr(var_count) ) // ' is' // &
    1381                                  'not available for radiation = .FALSE. or ' //&
    1382                                  'radiation_scheme = "constant"'
    1383                 CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
    1384              ELSE
    1385                 dopr_index(var_count) = 101
    1386                 dopr_unit  = 'W/m2'
    1387                 hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions+1 )
    1388                 unit = dopr_unit
    1389              ENDIF
    1390 
    1391           CASE ( 'rad_sw_in' )
    1392              IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
    1393              THEN
    1394                 message_string = 'data_output_pr = ' //                        &
    1395                                  TRIM( data_output_pr(var_count) ) // ' is' // &
    1396                                  'not available for radiation = .FALSE. or ' //&
    1397                                  'radiation_scheme = "constant"'
    1398                 CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
    1399              ELSE
    1400                 dopr_index(var_count) = 102
    1401                 dopr_unit  = 'W/m2'
    1402                 hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
    1403                 unit = dopr_unit
    1404              ENDIF
    1405 
    1406           CASE ( 'rad_sw_out')
    1407              IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
    1408              THEN
    1409                 message_string = 'data_output_pr = ' //                        &
    1410                                  TRIM( data_output_pr(var_count) ) // ' is' // &
    1411                                  'not available for radiation = .FALSE. or ' //&
    1412                                  'radiation_scheme = "constant"'
    1413                 CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
    1414              ELSE
    1415                 dopr_index(var_count) = 103
    1416                 dopr_unit  = 'W/m2'
    1417                 hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
    1418                 unit = dopr_unit
    1419              ENDIF
    1420 
    1421           CASE ( 'rad_lw_cs_hr' )
    1422              IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
    1423              THEN
    1424                 message_string = 'data_output_pr = ' //                        &
    1425                                  TRIM( data_output_pr(var_count) ) // ' is' // &
    1426                                  'not available for radiation = .FALSE. or ' //&
    1427                                  'radiation_scheme /= "rrtmg"'
    1428                 CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
    1429              ELSE
    1430                 dopr_index(var_count) = 104
    1431                 dopr_unit  = 'K/h'
    1432                 hom(:,2,104,:)  = SPREAD( zu, 2, statistic_regions+1 )
    1433                 unit = dopr_unit
    1434              ENDIF
    1435 
    1436           CASE ( 'rad_lw_hr' )
    1437              IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
    1438              THEN
    1439                 message_string = 'data_output_pr = ' //                        &
    1440                                  TRIM( data_output_pr(var_count) ) // ' is' // &
    1441                                  'not available for radiation = .FALSE. or ' //&
    1442                                  'radiation_scheme /= "rrtmg"'
    1443                 CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
    1444              ELSE
    1445                 dopr_index(var_count) = 105
    1446                 dopr_unit  = 'K/h'
    1447                 hom(:,2,105,:)  = SPREAD( zu, 2, statistic_regions+1 )
    1448                 unit = dopr_unit
    1449              ENDIF
    1450 
    1451           CASE ( 'rad_sw_cs_hr' )
    1452              IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
    1453              THEN
    1454                 message_string = 'data_output_pr = ' //                        &
    1455                                  TRIM( data_output_pr(var_count) ) // ' is' // &
    1456                                  'not available for radiation = .FALSE. or ' //&
    1457                                  'radiation_scheme /= "rrtmg"'
    1458                 CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
    1459              ELSE
    1460                 dopr_index(var_count) = 106
    1461                 dopr_unit  = 'K/h'
    1462                 hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions+1 )
    1463                 unit = dopr_unit
    1464              ENDIF
    1465 
    1466           CASE ( 'rad_sw_hr' )
    1467              IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
    1468              THEN
    1469                 message_string = 'data_output_pr = ' //                        &
    1470                                  TRIM( data_output_pr(var_count) ) // ' is' // &
    1471                                  'not available for radiation = .FALSE. or ' //&
    1472                                  'radiation_scheme /= "rrtmg"'
    1473                 CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
    1474              ELSE
    1475                 dopr_index(var_count) = 107
    1476                 dopr_unit  = 'K/h'
    1477                 hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions+1 )
    1478                 unit = dopr_unit
    1479              ENDIF
    1480 
    1481 
    1482           CASE DEFAULT
    1483              unit = 'illegal'
    1484 
    1485        END SELECT
    1486 
    1487 
    1488     END SUBROUTINE radiation_check_data_output_pr
    1489 
    1490 
    1491 !------------------------------------------------------------------------------!
     1519!--------------------------------------------------------------------------------------------------!
     1520 SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, dopr_unit )
     1521
     1522    USE arrays_3d,                                                                                 &
     1523        ONLY: zu
     1524
     1525    USE control_parameters,                                                                        &
     1526        ONLY: data_output_pr,                                                                      &
     1527              message_string
     1528
     1529    USE indices
     1530
     1531    USE profil_parameter
     1532
     1533    USE statistics
     1534
     1535    IMPLICIT NONE
     1536
     1537    CHARACTER(LEN=*) ::  unit       !<
     1538    CHARACTER(LEN=*) ::  variable   !<
     1539    CHARACTER(LEN=*) ::  dopr_unit  !< local value of dopr_unit
     1540
     1541    INTEGER(iwp) ::  var_count  !<
     1542
     1543    SELECT CASE ( TRIM( variable ) )
     1544
     1545      CASE ( 'rad_net' )
     1546          IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )  THEN
     1547             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
     1548                              'not available for radiation = .FALSE. or ' //                       &
     1549                              'radiation_scheme = "constant"'
     1550             CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
     1551          ELSE
     1552             dopr_index(var_count) = 99
     1553             dopr_unit  = 'W/m2'
     1554             hom(:,2,99,:)  = SPREAD( zw, 2, statistic_regions + 1 )
     1555             unit = dopr_unit
     1556          ENDIF
     1557
     1558       CASE ( 'rad_lw_in' )
     1559          IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' )  THEN
     1560             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
     1561                              'not available for radiation = .FALSE. or ' //                       &
     1562                              'radiation_scheme = "constant"'
     1563             CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
     1564          ELSE
     1565             dopr_index(var_count) = 100
     1566             dopr_unit  = 'W/m2'
     1567             hom(:,2,100,:)  = SPREAD( zw, 2, statistic_regions + 1 )
     1568             unit = dopr_unit
     1569          ENDIF
     1570
     1571       CASE ( 'rad_lw_out' )
     1572          IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' )  THEN
     1573             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
     1574                              'not available for radiation = .FALSE. or ' //                       &
     1575                              'radiation_scheme = "constant"'
     1576             CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
     1577          ELSE
     1578             dopr_index(var_count) = 101
     1579             dopr_unit  = 'W/m2'
     1580             hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions + 1 )
     1581             unit = dopr_unit
     1582          ENDIF
     1583
     1584       CASE ( 'rad_sw_in' )
     1585          IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' )  THEN
     1586             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
     1587                              'not available for radiation = .FALSE. or ' //                       &
     1588                              'radiation_scheme = "constant"'
     1589             CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
     1590          ELSE
     1591             dopr_index(var_count) = 102
     1592             dopr_unit  = 'W/m2'
     1593             hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions + 1 )
     1594             unit = dopr_unit
     1595          ENDIF
     1596
     1597       CASE ( 'rad_sw_out')
     1598          IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )  THEN
     1599             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
     1600                              'not available for radiation = .FALSE. or ' //                       &
     1601                              'radiation_scheme = "constant"'
     1602             CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
     1603          ELSE
     1604             dopr_index(var_count) = 103
     1605             dopr_unit  = 'W/m2'
     1606             hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions + 1 )
     1607             unit = dopr_unit
     1608          ENDIF
     1609
     1610       CASE ( 'rad_lw_cs_hr' )
     1611          IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )  THEN
     1612             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
     1613                              'not available for radiation = .FALSE. or ' //                       &
     1614                              'radiation_scheme /= "rrtmg"'
     1615             CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
     1616          ELSE
     1617             dopr_index(var_count) = 104
     1618             dopr_unit  = 'K/h'
     1619             hom(:,2,104,:)  = SPREAD( zu, 2, statistic_regions + 1 )
     1620             unit = dopr_unit
     1621          ENDIF
     1622
     1623       CASE ( 'rad_lw_hr' )
     1624          IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )  THEN
     1625             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
     1626                              'not available for radiation = .FALSE. or ' //                       &
     1627                              'radiation_scheme /= "rrtmg"'
     1628             CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
     1629          ELSE
     1630             dopr_index(var_count) = 105
     1631             dopr_unit  = 'K/h'
     1632             hom(:,2,105,:)  = SPREAD( zu, 2, statistic_regions + 1 )
     1633             unit = dopr_unit
     1634          ENDIF
     1635
     1636       CASE ( 'rad_sw_cs_hr' )
     1637          IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )  THEN
     1638             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
     1639                              'not available for radiation = .FALSE. or ' //                       &
     1640                              'radiation_scheme /= "rrtmg"'
     1641             CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
     1642          ELSE
     1643             dopr_index(var_count) = 106
     1644             dopr_unit  = 'K/h'
     1645             hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions + 1 )
     1646             unit = dopr_unit
     1647          ENDIF
     1648
     1649       CASE ( 'rad_sw_hr' )
     1650          IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )  THEN
     1651             message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // ' is' // &
     1652                              'not available for radiation = .FALSE. or ' //                       &
     1653                              'radiation_scheme /= "rrtmg"'
     1654             CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
     1655          ELSE
     1656             dopr_index(var_count) = 107
     1657             dopr_unit  = 'K/h'
     1658             hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions + 1 )
     1659             unit = dopr_unit
     1660          ENDIF
     1661
     1662
     1663       CASE DEFAULT
     1664          unit = 'illegal'
     1665
     1666    END SELECT
     1667
     1668
     1669 END SUBROUTINE radiation_check_data_output_pr
     1670
     1671
     1672!--------------------------------------------------------------------------------------------------!
    14921673! Description:
    14931674! ------------
    14941675!> Check parameters routine for radiation model
    1495 !------------------------------------------------------------------------------!
    1496     SUBROUTINE radiation_check_parameters
    1497 
    1498        USE control_parameters,                                                 &
    1499            ONLY: land_surface, message_string, urban_surface
    1500 
    1501        USE netcdf_data_input_mod,                                              &
    1502            ONLY:  input_pids_static
    1503 
    1504        IMPLICIT NONE
    1505 
    1506 !
    1507 !--    In case no urban-surface or land-surface model is applied, usage of
    1508 !--    a radiation model make no sense.
    1509        IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
    1510           message_string = 'Usage of radiation module is only allowed if ' //  &
    1511                            'land-surface and/or urban-surface model is applied.'
    1512           CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
    1513        ENDIF
    1514 
    1515        IF ( radiation_scheme /= 'constant'   .AND.                             &
    1516             radiation_scheme /= 'clear-sky'  .AND.                             &
    1517             radiation_scheme /= 'rrtmg'      .AND.                             &
    1518             radiation_scheme /= 'external' )  THEN
    1519           message_string = 'unknown radiation_scheme = '//                     &
    1520                            TRIM( radiation_scheme )
    1521           CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
    1522        ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
     1676!--------------------------------------------------------------------------------------------------!
     1677 SUBROUTINE radiation_check_parameters
     1678
     1679    USE control_parameters,                                                                        &
     1680        ONLY: land_surface,                                                                        &
     1681              message_string,                                                                      &
     1682              urban_surface
     1683
     1684    USE netcdf_data_input_mod,                                                                     &
     1685        ONLY:  input_pids_static
     1686
     1687    IMPLICIT NONE
     1688
     1689!
     1690!-- In case no urban-surface or land-surface model is applied, usage of a radiation model makes no
     1691!-- sense.
     1692    IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
     1693       message_string = 'Usage of radiation module is only allowed if ' //                         &
     1694                        'land-surface and/or urban-surface model is applied.'
     1695       CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
     1696    ENDIF
     1697
     1698    IF ( radiation_scheme /= 'constant'  .AND.  radiation_scheme /= 'clear-sky'  .AND.             &
     1699         radiation_scheme /= 'rrtmg'  .AND.  radiation_scheme /= 'external' )  THEN
     1700       message_string = 'unknown radiation_scheme = '// TRIM( radiation_scheme )
     1701       CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
     1702    ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
    15231703#if ! defined ( __rrtmg )
    1524           message_string = 'radiation_scheme = "rrtmg" requires ' //           &
    1525                            'compilation of PALM with pre-processor ' //        &
    1526                            'directive -D__rrtmg'
    1527           CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
     1704       message_string = 'radiation_scheme = "rrtmg" requires compilation of PALM with ' //         &
     1705                        'pre-processor directive -D__rrtmg'
     1706       CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
    15281707#endif
    15291708#if defined ( __rrtmg ) && ! defined( __netcdf )
    1530           message_string = 'radiation_scheme = "rrtmg" requires ' //           &
    1531                            'the use of NetCDF (preprocessor directive ' //     &
    1532                            '-D__netcdf'
    1533           CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
     1709       message_string = 'radiation_scheme = "rrtmg" requires the use of NetCDF ' //                &
     1710                        '(preprocessor directive -D__netcdf'
     1711       CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
    15341712#endif
    15351713
    1536        ENDIF
    1537 !
    1538 !--    Checks performed only if data is given via namelist only.
    1539        IF ( .NOT. input_pids_static )  THEN
    1540           IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.          &
    1541                radiation_scheme == 'clear-sky')  THEN
    1542              message_string = 'radiation_scheme = "clear-sky" in combination'//&
    1543                               'with albedo_type = 0 requires setting of'//     &
    1544                               'albedo /= 9999999.9'
    1545              CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
    1546           ENDIF
    1547 
    1548           IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.     &
    1549              ( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
    1550           .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp&
    1551              ) ) THEN
    1552              message_string = 'radiation_scheme = "rrtmg" in combination' //   &
    1553                               'with albedo_type = 0 requires setting of ' //   &
    1554                               'albedo_lw_dif /= 9999999.9' //                  &
    1555                               'albedo_lw_dir /= 9999999.9' //                  &
    1556                               'albedo_sw_dif /= 9999999.9 and' //              &
    1557                               'albedo_sw_dir /= 9999999.9'
    1558              CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
    1559           ENDIF
    1560        ENDIF
    1561 !
    1562 !--    Parallel rad_angular_discretization without raytrace_mpi_rma is not implemented
    1563 !--    Serial mode does not allow mpi_rma
     1714    ENDIF
     1715!
     1716!-- Checks performed if data is given via namelist only.
     1717    IF ( .NOT. input_pids_static )  THEN
     1718       IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.                                 &
     1719            radiation_scheme == 'clear-sky')  THEN
     1720          message_string = 'radiation_scheme = "clear-sky" in combination with albedo_type = 0 ' //&
     1721                           'requires setting of albedo /= 9999999.9'
     1722          CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
     1723       ENDIF
     1724
     1725       IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.                            &
     1726          ( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp                       &
     1727       .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp ) )  THEN
     1728          message_string = 'radiation_scheme = "rrtmg" in combination with albedo_type = 0 ' //    &
     1729                           'requires setting of albedo_lw_dif /= 9999999.9' //                     &
     1730                           'albedo_lw_dir /= 9999999.9' //                                         &
     1731                           'albedo_sw_dif /= 9999999.9 and albedo_sw_dir /= 9999999.9'
     1732          CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
     1733       ENDIF
     1734    ENDIF
     1735!
     1736!-- Parallel rad_angular_discretization without raytrace_mpi_rma is not implemented. Serial mode does
     1737!-- not allow mpi_rma
    15641738#if defined( __parallel )
    1565        IF ( rad_angular_discretization  .AND.  .NOT. raytrace_mpi_rma )  THEN
    1566           message_string = 'rad_angular_discretization can only be used ' //  &
    1567                            'together with raytrace_mpi_rma or when ' //  &
    1568                            'no parallelization is applied.'
    1569           CALL message( 'readiation_check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
    1570        ENDIF
     1739    IF ( rad_angular_discretization  .AND.  .NOT. raytrace_mpi_rma )  THEN
     1740       message_string = 'rad_angular_discretization can only be used together with ' //            &
     1741                        'raytrace_mpi_rma or when no parallelization is applied.'
     1742       CALL message( 'readiation_check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
     1743    ENDIF
    15711744#else
    1572        IF ( raytrace_mpi_rma )  THEN
    1573           message_string = 'raytrace_mpi_rma = .T. not allowed in serial mode'
    1574           CALL message( 'readiation_check_parameters', 'PA0710', 1, 2, 0, 6, 0 )
    1575        ENDIF
     1745    IF ( raytrace_mpi_rma )  THEN
     1746       message_string = 'raytrace_mpi_rma = .T. not allowed in serial mode'
     1747       CALL message( 'readiation_check_parameters', 'PA0710', 1, 2, 0, 6, 0 )
     1748    ENDIF
    15761749#endif
    15771750
    1578        IF ( cloud_droplets  .AND.   radiation_scheme == 'rrtmg'  .AND.         &
    1579             average_radiation ) THEN
    1580           message_string = 'average_radiation = .T. with radiation_scheme'//   &
    1581                            '= "rrtmg" in combination cloud_droplets = .T.'//   &
    1582                            'is not implementd'
    1583           CALL message( 'check_parameters', 'PA0560', 1, 2, 0, 6, 0 )
    1584        ENDIF
    1585 
    1586 !
    1587 !--    Incialize svf normalization reporting histogram
    1588        svfnorm_report_num = 1
    1589        DO WHILE ( svfnorm_report_thresh(svfnorm_report_num) < 1e20_wp          &
    1590                    .AND. svfnorm_report_num <= 30 )
    1591           svfnorm_report_num = svfnorm_report_num + 1
    1592        ENDDO
    1593        svfnorm_report_num = svfnorm_report_num - 1
    1594 !
    1595 !--    Check for dt_radiation
    1596        IF ( dt_radiation <= 0.0 )  THEN
    1597           message_string = 'dt_radiation must be > 0.0'
    1598           CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 )
    1599        ENDIF
    1600 
    1601     END SUBROUTINE radiation_check_parameters
    1602 
    1603 
    1604 !------------------------------------------------------------------------------!
     1751    IF ( cloud_droplets  .AND.   radiation_scheme == 'rrtmg'  .AND.  average_radiation )  THEN
     1752       message_string = 'average_radiation = .T. with radiation_scheme = "rrtmg" ' //              &
     1753                        'in combination cloud_droplets = .T. is not implementd'
     1754       CALL message( 'check_parameters', 'PA0560', 1, 2, 0, 6, 0 )
     1755    ENDIF
     1756
     1757!
     1758!-- Initialize svf normalization reporting histogram
     1759    svfnorm_report_num = 1
     1760    DO WHILE ( svfnorm_report_thresh(svfnorm_report_num) < 1E20_wp .AND. svfnorm_report_num <= 30 )
     1761       svfnorm_report_num = svfnorm_report_num + 1
     1762    ENDDO
     1763    svfnorm_report_num = svfnorm_report_num - 1
     1764!
     1765!-- Check for dt_radiation
     1766    IF ( dt_radiation <= 0.0 )  THEN
     1767       message_string = 'dt_radiation must be > 0.0'
     1768       CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 )
     1769    ENDIF
     1770
     1771 END SUBROUTINE radiation_check_parameters
     1772
     1773
     1774!--------------------------------------------------------------------------------------------------!
    16051775! Description:
    16061776! ------------
    16071777!> Initialization of the radiation model and Radiative Transfer Model
    1608 !------------------------------------------------------------------------------!
    1609     SUBROUTINE radiation_init
    1610 
    1611        IMPLICIT NONE
    1612 
    1613        INTEGER(iwp) ::  i         !< running index x-direction
    1614        INTEGER(iwp) ::  is        !< running index for input surface elements
    1615        INTEGER(iwp) ::  ioff      !< offset in x between surface element reference grid point in atmosphere and actual surface
    1616        INTEGER(iwp) ::  j         !< running index y-direction
    1617        INTEGER(iwp) ::  joff      !< offset in y between surface element reference grid point in atmosphere and actual surface
    1618        INTEGER(iwp) ::  k         !< running index z-direction
    1619        INTEGER(iwp) ::  l         !< running index for orientation of vertical surfaces
    1620        INTEGER(iwp) ::  m         !< running index for surface elements
    1621        INTEGER(iwp) ::  ntime = 0 !< number of available external radiation timesteps
     1778!--------------------------------------------------------------------------------------------------!
     1779 SUBROUTINE radiation_init
     1780
     1781    IMPLICIT NONE
     1782
     1783    INTEGER(iwp) ::  i          !< running index x-direction
     1784    INTEGER(iwp) ::  is         !< running index for input surface elements
     1785    INTEGER(iwp) ::  ioff       !< offset in x between surface element reference grid point in atmosphere and actual surface
     1786    INTEGER(iwp) ::  j          !< running index y-direction
     1787    INTEGER(iwp) ::  joff       !< offset in y between surface element reference grid point in atmosphere and actual surface
     1788    INTEGER(iwp) ::  k          !< running index z-direction
     1789    INTEGER(iwp) ::  l          !< running index for orientation of vertical surfaces
     1790    INTEGER(iwp) ::  m          !< running index for surface elements
     1791    INTEGER(iwp) ::  ntime = 0 !< number of available external radiation timesteps
    16221792#if defined( __rrtmg )
    1623        INTEGER(iwp) ::  ind_type  !< running index for subgrid-surface tiles
     1793    INTEGER(iwp) :: ind_type   !< running index for subgrid-surface tiles
    16241794#endif
    1625        LOGICAL      ::  radiation_input_root_domain !< flag indicating the existence of a dynamic input file for the root domain
    1626 
    1627 
    1628        IF ( debug_output )  CALL debug_message( 'radiation_init', 'start' )
    1629 !
    1630 !--    Activate radiation_interactions according to the existence of vertical surfaces and/or trees
    1631 !      or if biometeorology output is required for flat surfaces.
    1632 !--    The namelist parameter radiation_interactions_on can override this behavior.
    1633 !--    (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in
    1634 !--    init_surface_arrays.)
    1635        IF ( radiation_interactions_on )  THEN
    1636           IF ( vertical_surfaces_exist  .OR.  plant_canopy  .OR.  biometeorology )  THEN
    1637              radiation_interactions    = .TRUE.
    1638              average_radiation         = .TRUE.
    1639           ELSE
    1640              radiation_interactions_on = .FALSE.   !< reset namelist parameter: no interactions
    1641                                                    !< calculations necessary in case of flat surface
     1795    LOGICAL ::  radiation_input_root_domain  !< flag indicating the existence of a dynamic input file for the root domain
     1796
     1797
     1798    IF ( debug_output )  CALL debug_message( 'radiation_init', 'start' )
     1799!
     1800!-- Activate radiation_interactions according to the existence of vertical surfaces and/or trees
     1801!   or if biometeorology output is required for flat surfaces.
     1802!-- The namelist parameter radiation_interactions_on can override this behavior. (This check cannot
     1803!-- be performed in check_parameters, because vertical_surfaces_exist is first set in
     1804!-- init_surface_arrays.)
     1805    IF ( radiation_interactions_on )  THEN
     1806       IF ( vertical_surfaces_exist  .OR.  plant_canopy  .OR.  biometeorology )  THEN
     1807          radiation_interactions    = .TRUE.
     1808          average_radiation         = .TRUE.
     1809       ELSE
     1810          radiation_interactions_on = .FALSE.  !< reset namelist parameter: no interactions
     1811                                               !< calculations necessary in case of flat surface
     1812       ENDIF
     1813    ELSEIF ( vertical_surfaces_exist  .OR.  plant_canopy  .OR.  biometeorology )  THEN
     1814       message_string = 'radiation_interactions_on is set to .FALSE. although vertical ' //        &
     1815                        'surfaces and/or trees or biometeorology exist is ON. The model will ' //  &
     1816                        'run without RTM (no shadows, no radiation reflections)'
     1817       CALL message( 'init_3d_model', 'PA0348', 0, 1, 0, 6, 0 )
     1818    ENDIF
     1819!
     1820!-- Precalculate some time constants
     1821    d_hours_day    = 1.0_wp / REAL( hours_per_day, KIND = wp )
     1822    d_seconds_hour = 1.0_wp / seconds_per_hour
     1823
     1824!
     1825!-- If required, initialize radiation interactions between surfaces via sky-view factors. This must
     1826!-- be done before radiation is initialized.
     1827    IF ( radiation_interactions )  CALL radiation_interaction_init
     1828!
     1829!-- Allocate array for storing the surface net radiation
     1830    IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_net )  .AND.  surf_lsm_h%ns > 0 )  THEN
     1831       ALLOCATE( surf_lsm_h%rad_net(1:surf_lsm_h%ns) )
     1832       surf_lsm_h%rad_net = 0.0_wp
     1833    ENDIF
     1834    IF ( .NOT. ALLOCATED ( surf_usm_h%rad_net )  .AND.  surf_usm_h%ns > 0 )  THEN
     1835       ALLOCATE( surf_usm_h%rad_net(1:surf_usm_h%ns) )
     1836       surf_usm_h%rad_net = 0.0_wp
     1837    ENDIF
     1838    DO  l = 0, 3
     1839       IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_net )  .AND.  surf_lsm_v(l)%ns > 0 )  THEN
     1840          ALLOCATE( surf_lsm_v(l)%rad_net(1:surf_lsm_v(l)%ns) )
     1841          surf_lsm_v(l)%rad_net = 0.0_wp
     1842       ENDIF
     1843       IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_net )  .AND.  surf_usm_v(l)%ns > 0 )  THEN
     1844          ALLOCATE( surf_usm_v(l)%rad_net(1:surf_usm_v(l)%ns) )
     1845          surf_usm_v(l)%rad_net = 0.0_wp
     1846       ENDIF
     1847    ENDDO
     1848
     1849
     1850!
     1851!-- Allocate array for storing the surface longwave (out) radiation change
     1852    IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_lw_out_change_0 )  .AND.  surf_lsm_h%ns > 0 )  THEN
     1853       ALLOCATE( surf_lsm_h%rad_lw_out_change_0(1:surf_lsm_h%ns) )
     1854       surf_lsm_h%rad_lw_out_change_0 = 0.0_wp
     1855    ENDIF
     1856    IF ( .NOT. ALLOCATED ( surf_usm_h%rad_lw_out_change_0 )  .AND.  surf_usm_h%ns > 0 ) THEN
     1857       ALLOCATE( surf_usm_h%rad_lw_out_change_0(1:surf_usm_h%ns) )
     1858       surf_usm_h%rad_lw_out_change_0 = 0.0_wp
     1859    ENDIF
     1860    DO  l = 0, 3
     1861       IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_lw_out_change_0 )  .AND.  surf_lsm_v(l)%ns > 0 )   &
     1862       THEN
     1863          ALLOCATE( surf_lsm_v(l)%rad_lw_out_change_0(1:surf_lsm_v(l)%ns) )
     1864          surf_lsm_v(l)%rad_lw_out_change_0 = 0.0_wp
     1865       ENDIF
     1866       IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_lw_out_change_0 )  .AND.  surf_usm_v(l)%ns > 0 )   &
     1867       THEN
     1868          ALLOCATE( surf_usm_v(l)%rad_lw_out_change_0(1:surf_usm_v(l)%ns) )
     1869          surf_usm_v(l)%rad_lw_out_change_0 = 0.0_wp
     1870       ENDIF
     1871    ENDDO
     1872
     1873!
     1874!-- Allocate surface arrays for incoming/outgoing short/longwave radiation
     1875    IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_sw_in )  .AND.  surf_lsm_h%ns > 0 )  THEN
     1876       ALLOCATE( surf_lsm_h%rad_sw_in(1:surf_lsm_h%ns)  )
     1877       ALLOCATE( surf_lsm_h%rad_sw_out(1:surf_lsm_h%ns) )
     1878       ALLOCATE( surf_lsm_h%rad_sw_dir(1:surf_lsm_h%ns) )
     1879       ALLOCATE( surf_lsm_h%rad_sw_dif(1:surf_lsm_h%ns) )
     1880       ALLOCATE( surf_lsm_h%rad_sw_ref(1:surf_lsm_h%ns) )
     1881       ALLOCATE( surf_lsm_h%rad_sw_res(1:surf_lsm_h%ns) )
     1882       ALLOCATE( surf_lsm_h%rad_lw_in(1:surf_lsm_h%ns)  )
     1883       ALLOCATE( surf_lsm_h%rad_lw_out(1:surf_lsm_h%ns) )
     1884       ALLOCATE( surf_lsm_h%rad_lw_dif(1:surf_lsm_h%ns) )
     1885       ALLOCATE( surf_lsm_h%rad_lw_ref(1:surf_lsm_h%ns) )
     1886       ALLOCATE( surf_lsm_h%rad_lw_res(1:surf_lsm_h%ns) )
     1887       surf_lsm_h%rad_sw_in  = 0.0_wp
     1888       surf_lsm_h%rad_sw_out = 0.0_wp
     1889       surf_lsm_h%rad_sw_dir = 0.0_wp
     1890       surf_lsm_h%rad_sw_dif = 0.0_wp
     1891       surf_lsm_h%rad_sw_ref = 0.0_wp
     1892       surf_lsm_h%rad_sw_res = 0.0_wp
     1893       surf_lsm_h%rad_lw_in  = 0.0_wp
     1894       surf_lsm_h%rad_lw_out = 0.0_wp
     1895       surf_lsm_h%rad_lw_dif = 0.0_wp
     1896       surf_lsm_h%rad_lw_ref = 0.0_wp
     1897       surf_lsm_h%rad_lw_res = 0.0_wp
     1898    ENDIF
     1899    IF ( .NOT. ALLOCATED ( surf_usm_h%rad_sw_in )  .AND.  surf_usm_h%ns > 0 ) THEN
     1900       ALLOCATE( surf_usm_h%rad_sw_in(1:surf_usm_h%ns)  )
     1901       ALLOCATE( surf_usm_h%rad_sw_out(1:surf_usm_h%ns) )
     1902       ALLOCATE( surf_usm_h%rad_sw_dir(1:surf_usm_h%ns) )
     1903       ALLOCATE( surf_usm_h%rad_sw_dif(1:surf_usm_h%ns) )
     1904       ALLOCATE( surf_usm_h%rad_sw_ref(1:surf_usm_h%ns) )
     1905       ALLOCATE( surf_usm_h%rad_sw_res(1:surf_usm_h%ns) )
     1906       ALLOCATE( surf_usm_h%rad_lw_in(1:surf_usm_h%ns)  )
     1907       ALLOCATE( surf_usm_h%rad_lw_out(1:surf_usm_h%ns) )
     1908       ALLOCATE( surf_usm_h%rad_lw_dif(1:surf_usm_h%ns) )
     1909       ALLOCATE( surf_usm_h%rad_lw_ref(1:surf_usm_h%ns) )
     1910       ALLOCATE( surf_usm_h%rad_lw_res(1:surf_usm_h%ns) )
     1911       surf_usm_h%rad_sw_in  = 0.0_wp
     1912       surf_usm_h%rad_sw_out = 0.0_wp
     1913       surf_usm_h%rad_sw_dir = 0.0_wp
     1914       surf_usm_h%rad_sw_dif = 0.0_wp
     1915       surf_usm_h%rad_sw_ref = 0.0_wp
     1916       surf_usm_h%rad_sw_res = 0.0_wp
     1917       surf_usm_h%rad_lw_in  = 0.0_wp
     1918       surf_usm_h%rad_lw_out = 0.0_wp
     1919       surf_usm_h%rad_lw_dif = 0.0_wp
     1920       surf_usm_h%rad_lw_ref = 0.0_wp
     1921       surf_usm_h%rad_lw_res = 0.0_wp
     1922    ENDIF
     1923    DO  l = 0, 3
     1924       IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_sw_in )  .AND.  surf_lsm_v(l)%ns > 0 )  THEN
     1925          ALLOCATE( surf_lsm_v(l)%rad_sw_in(1:surf_lsm_v(l)%ns)  )
     1926          ALLOCATE( surf_lsm_v(l)%rad_sw_out(1:surf_lsm_v(l)%ns) )
     1927          ALLOCATE( surf_lsm_v(l)%rad_sw_dir(1:surf_lsm_v(l)%ns) )
     1928          ALLOCATE( surf_lsm_v(l)%rad_sw_dif(1:surf_lsm_v(l)%ns) )
     1929          ALLOCATE( surf_lsm_v(l)%rad_sw_ref(1:surf_lsm_v(l)%ns) )
     1930          ALLOCATE( surf_lsm_v(l)%rad_sw_res(1:surf_lsm_v(l)%ns) )
     1931
     1932          ALLOCATE( surf_lsm_v(l)%rad_lw_in(1:surf_lsm_v(l)%ns)  )
     1933          ALLOCATE( surf_lsm_v(l)%rad_lw_out(1:surf_lsm_v(l)%ns) )
     1934          ALLOCATE( surf_lsm_v(l)%rad_lw_dif(1:surf_lsm_v(l)%ns) )
     1935          ALLOCATE( surf_lsm_v(l)%rad_lw_ref(1:surf_lsm_v(l)%ns) )
     1936          ALLOCATE( surf_lsm_v(l)%rad_lw_res(1:surf_lsm_v(l)%ns) )
     1937
     1938          surf_lsm_v(l)%rad_sw_in  = 0.0_wp
     1939          surf_lsm_v(l)%rad_sw_out = 0.0_wp
     1940          surf_lsm_v(l)%rad_sw_dir = 0.0_wp
     1941          surf_lsm_v(l)%rad_sw_dif = 0.0_wp
     1942          surf_lsm_v(l)%rad_sw_ref = 0.0_wp
     1943          surf_lsm_v(l)%rad_sw_res = 0.0_wp
     1944
     1945          surf_lsm_v(l)%rad_lw_in  = 0.0_wp
     1946          surf_lsm_v(l)%rad_lw_out = 0.0_wp
     1947          surf_lsm_v(l)%rad_lw_dif = 0.0_wp
     1948          surf_lsm_v(l)%rad_lw_ref = 0.0_wp
     1949          surf_lsm_v(l)%rad_lw_res = 0.0_wp
     1950       ENDIF
     1951       IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_sw_in )  .AND.  surf_usm_v(l)%ns > 0 )  THEN
     1952          ALLOCATE( surf_usm_v(l)%rad_sw_in(1:surf_usm_v(l)%ns)  )
     1953          ALLOCATE( surf_usm_v(l)%rad_sw_out(1:surf_usm_v(l)%ns) )
     1954          ALLOCATE( surf_usm_v(l)%rad_sw_dir(1:surf_usm_v(l)%ns) )
     1955          ALLOCATE( surf_usm_v(l)%rad_sw_dif(1:surf_usm_v(l)%ns) )
     1956          ALLOCATE( surf_usm_v(l)%rad_sw_ref(1:surf_usm_v(l)%ns) )
     1957          ALLOCATE( surf_usm_v(l)%rad_sw_res(1:surf_usm_v(l)%ns) )
     1958          ALLOCATE( surf_usm_v(l)%rad_lw_in(1:surf_usm_v(l)%ns)  )
     1959          ALLOCATE( surf_usm_v(l)%rad_lw_out(1:surf_usm_v(l)%ns) )
     1960          ALLOCATE( surf_usm_v(l)%rad_lw_dif(1:surf_usm_v(l)%ns) )
     1961          ALLOCATE( surf_usm_v(l)%rad_lw_ref(1:surf_usm_v(l)%ns) )
     1962          ALLOCATE( surf_usm_v(l)%rad_lw_res(1:surf_usm_v(l)%ns) )
     1963          surf_usm_v(l)%rad_sw_in  = 0.0_wp
     1964          surf_usm_v(l)%rad_sw_out = 0.0_wp
     1965          surf_usm_v(l)%rad_sw_dir = 0.0_wp
     1966          surf_usm_v(l)%rad_sw_dif = 0.0_wp
     1967          surf_usm_v(l)%rad_sw_ref = 0.0_wp
     1968          surf_usm_v(l)%rad_sw_res = 0.0_wp
     1969          surf_usm_v(l)%rad_lw_in  = 0.0_wp
     1970          surf_usm_v(l)%rad_lw_out = 0.0_wp
     1971          surf_usm_v(l)%rad_lw_dif = 0.0_wp
     1972          surf_usm_v(l)%rad_lw_ref = 0.0_wp
     1973          surf_usm_v(l)%rad_lw_res = 0.0_wp
     1974       ENDIF
     1975    ENDDO
     1976!
     1977!-- Fix net radiation in case of radiation_scheme = 'constant'
     1978    IF ( radiation_scheme == 'constant' )  THEN
     1979       IF ( ALLOCATED( surf_lsm_h%rad_net ) )  surf_lsm_h%rad_net    = net_radiation
     1980       IF ( ALLOCATED( surf_usm_h%rad_net ) )  surf_usm_h%rad_net    = net_radiation
     1981!
     1982!--    Todo: weight with inclination angle
     1983       DO  l = 0, 3
     1984          IF ( ALLOCATED( surf_lsm_v(l)%rad_net ) )  surf_lsm_v(l)%rad_net = net_radiation
     1985          IF ( ALLOCATED( surf_usm_v(l)%rad_net ) )  surf_usm_v(l)%rad_net = net_radiation
     1986       ENDDO
     1987!       radiation = .FALSE.
     1988!
     1989!-- Calculate orbital constants
     1990    ELSE
     1991       decl_1 = SIN( 23.45_wp * pi / 180.0_wp )
     1992       decl_2 = 2.0_wp * pi / 365.0_wp
     1993       decl_3 = decl_2 * 81.0_wp
     1994       lat    = latitude * pi / 180.0_wp
     1995       lon    = longitude * pi / 180.0_wp
     1996    ENDIF
     1997
     1998    IF ( radiation_scheme == 'clear-sky'  .OR.                                                     &
     1999         radiation_scheme == 'constant'   .OR.                                                     &
     2000         radiation_scheme == 'external' )  THEN
     2001!
     2002!--    Allocate arrays for incoming/outgoing short/longwave radiation
     2003       IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
     2004          ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
     2005          rad_sw_in = 0.0_wp
     2006       ENDIF
     2007       IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
     2008          ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
     2009          rad_sw_out = 0.0_wp
     2010       ENDIF
     2011
     2012       IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
     2013          ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
     2014          rad_lw_in = 0.0_wp
     2015       ENDIF
     2016       IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
     2017          ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
     2018          rad_lw_out = 0.0_wp
     2019       ENDIF
     2020
     2021!
     2022!--    Allocate average arrays for incoming/outgoing short/longwave radiation
     2023       IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
     2024          ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
     2025       ENDIF
     2026       IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
     2027          ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
     2028       ENDIF
     2029
     2030       IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
     2031          ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
     2032       ENDIF
     2033       IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
     2034          ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
     2035       ENDIF
     2036!
     2037!--    Allocate arrays for broadband albedo, and level 1 initialization via namelist paramter,
     2038!--    unless not already allocated.
     2039       IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )  THEN
     2040          ALLOCATE( surf_lsm_h%albedo(1:surf_lsm_h%ns,0:2) )
     2041          surf_lsm_h%albedo    = albedo
     2042       ENDIF
     2043       IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )  THEN
     2044          ALLOCATE( surf_usm_h%albedo(1:surf_usm_h%ns,0:2) )
     2045          surf_usm_h%albedo    = albedo
     2046       ENDIF
     2047
     2048       DO  l = 0, 3
     2049          IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )  THEN
     2050             ALLOCATE( surf_lsm_v(l)%albedo(1:surf_lsm_v(l)%ns,0:2) )
     2051             surf_lsm_v(l)%albedo = albedo
    16422052          ENDIF
    1643        ELSEIF ( vertical_surfaces_exist  .OR.  plant_canopy  .OR.  biometeorology )  THEN
    1644           message_string = 'radiation_interactions_on is set to .FALSE. although '     // &
    1645                            'vertical surfaces and/or trees or biometeorology exist '   // &
    1646                            'is ON. The model will run without RTM (no shadows, no '    // &
    1647                            'radiation reflections)'
    1648           CALL message( 'init_3d_model', 'PA0348', 0, 1, 0, 6, 0 )
    1649        ENDIF
    1650 !
    1651 !--    Precalculate some time constants
    1652        d_hours_day    = 1.0_wp / REAL( hours_per_day, KIND = wp )
    1653        d_seconds_hour = 1.0_wp / seconds_per_hour
    1654 
    1655 !
    1656 !--    If required, initialize radiation interactions between surfaces
    1657 !--    via sky-view factors. This must be done before radiation is initialized.
    1658        IF ( radiation_interactions )  CALL radiation_interaction_init
    1659 !
    1660 !--    Allocate array for storing the surface net radiation
    1661        IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_net )  .AND.                      &
    1662                   surf_lsm_h%ns > 0  )   THEN
    1663           ALLOCATE( surf_lsm_h%rad_net(1:surf_lsm_h%ns) )
    1664           surf_lsm_h%rad_net = 0.0_wp
    1665        ENDIF
    1666        IF ( .NOT. ALLOCATED ( surf_usm_h%rad_net )  .AND.                      &
    1667                   surf_usm_h%ns > 0  )  THEN
    1668           ALLOCATE( surf_usm_h%rad_net(1:surf_usm_h%ns) )
    1669           surf_usm_h%rad_net = 0.0_wp
    1670        ENDIF
    1671        DO  l = 0, 3
    1672           IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_net )  .AND.                &
    1673                      surf_lsm_v(l)%ns > 0  )  THEN
    1674              ALLOCATE( surf_lsm_v(l)%rad_net(1:surf_lsm_v(l)%ns) )
    1675              surf_lsm_v(l)%rad_net = 0.0_wp
    1676           ENDIF
    1677           IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_net )  .AND.                &
    1678                      surf_usm_v(l)%ns > 0  )  THEN
    1679              ALLOCATE( surf_usm_v(l)%rad_net(1:surf_usm_v(l)%ns) )
    1680              surf_usm_v(l)%rad_net = 0.0_wp
     2053          IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )  THEN
     2054             ALLOCATE( surf_usm_v(l)%albedo(1:surf_usm_v(l)%ns,0:2) )
     2055             surf_usm_v(l)%albedo = albedo
    16812056          ENDIF
    16822057       ENDDO
    1683 
    1684 
    1685 !
    1686 !--    Allocate array for storing the surface longwave (out) radiation change
    1687        IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_lw_out_change_0 )  .AND.          &
    1688                   surf_lsm_h%ns > 0  )   THEN
    1689           ALLOCATE( surf_lsm_h%rad_lw_out_change_0(1:surf_lsm_h%ns) )
    1690           surf_lsm_h%rad_lw_out_change_0 = 0.0_wp
    1691        ENDIF
    1692        IF ( .NOT. ALLOCATED ( surf_usm_h%rad_lw_out_change_0 )  .AND.          &
    1693                   surf_usm_h%ns > 0  )  THEN
    1694           ALLOCATE( surf_usm_h%rad_lw_out_change_0(1:surf_usm_h%ns) )
    1695           surf_usm_h%rad_lw_out_change_0 = 0.0_wp
    1696        ENDIF
     2058!
     2059!--    Level 2 initialization of broadband albedo via given albedo_type.
     2060!--    Only if albedo_type is non-zero. In case of urban surface and input data is read from ASCII
     2061!--    file, albedo_type will be zero, so that albedo won't be overwritten.
     2062       DO  m = 1, surf_lsm_h%ns
     2063          IF ( surf_lsm_h%albedo_type(m,ind_veg_wall) /= 0 )                                       &
     2064             surf_lsm_h%albedo(m,ind_veg_wall) =                                                   &
     2065             albedo_pars(0,surf_lsm_h%albedo_type(m,ind_veg_wall))
     2066          IF ( surf_lsm_h%albedo_type(m,ind_pav_green) /= 0 )                                      &
     2067             surf_lsm_h%albedo(m,ind_pav_green) =                                                  &
     2068                        albedo_pars(0,surf_lsm_h%albedo_type(m,ind_pav_green))
     2069          IF ( surf_lsm_h%albedo_type(m,ind_wat_win) /= 0 )                                        &
     2070             surf_lsm_h%albedo(m,ind_wat_win) =                                                    &
     2071                        albedo_pars(0,surf_lsm_h%albedo_type(m,ind_wat_win))
     2072       ENDDO
     2073       DO  m = 1, surf_usm_h%ns
     2074          IF ( surf_usm_h%albedo_type(m,ind_veg_wall) /= 0 )                                       &
     2075             surf_usm_h%albedo(m,ind_veg_wall) =                                                   &
     2076                        albedo_pars(0,surf_usm_h%albedo_type(m,ind_veg_wall))
     2077          IF ( surf_usm_h%albedo_type(m,ind_pav_green) /= 0 )                                      &
     2078             surf_usm_h%albedo(m,ind_pav_green) =                                                  &
     2079                        albedo_pars(0,surf_usm_h%albedo_type(m,ind_pav_green))
     2080          IF ( surf_usm_h%albedo_type(m,ind_wat_win) /= 0 )                                        &
     2081             surf_usm_h%albedo(m,ind_wat_win) =                                                    &
     2082                        albedo_pars(0,surf_usm_h%albedo_type(m,ind_wat_win))
     2083       ENDDO
     2084
    16972085       DO  l = 0, 3
    1698           IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_lw_out_change_0 )  .AND.    &
    1699                      surf_lsm_v(l)%ns > 0  )  THEN
    1700              ALLOCATE( surf_lsm_v(l)%rad_lw_out_change_0(1:surf_lsm_v(l)%ns) )
    1701              surf_lsm_v(l)%rad_lw_out_change_0 = 0.0_wp
     2086          DO  m = 1, surf_lsm_v(l)%ns
     2087             IF ( surf_lsm_v(l)%albedo_type(m,ind_veg_wall) /= 0 )                                 &
     2088                surf_lsm_v(l)%albedo(m,ind_veg_wall) =                                             &
     2089                     albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_veg_wall))
     2090             IF ( surf_lsm_v(l)%albedo_type(m,ind_pav_green) /= 0 )                                &
     2091                surf_lsm_v(l)%albedo(m,ind_pav_green) =                                            &
     2092                     albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_pav_green))
     2093             IF ( surf_lsm_v(l)%albedo_type(m,ind_wat_win) /= 0 )                                  &
     2094                surf_lsm_v(l)%albedo(m,ind_wat_win) =                                              &
     2095                     albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_wat_win))
     2096          ENDDO
     2097          DO  m = 1, surf_usm_v(l)%ns
     2098             IF ( surf_usm_v(l)%albedo_type(m,ind_veg_wall) /= 0 )                                 &
     2099                surf_usm_v(l)%albedo(m,ind_veg_wall) =                                             &
     2100                     albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_veg_wall))
     2101             IF ( surf_usm_v(l)%albedo_type(m,ind_pav_green) /= 0 )                                &
     2102                surf_usm_v(l)%albedo(m,ind_pav_green) =                                            &
     2103                     albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_pav_green))
     2104             IF ( surf_usm_v(l)%albedo_type(m,ind_wat_win) /= 0 )                                  &
     2105                surf_usm_v(l)%albedo(m,ind_wat_win) =                                              &
     2106                     albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_wat_win))
     2107          ENDDO
     2108       ENDDO
     2109
     2110!
     2111!--    Level 3 initialization at grid points where albedo type is zero.
     2112!--    This case, albedo is taken from file. In case of constant radiation or clear sky, only
     2113!--    broadband albedo is given.
     2114       IF ( albedo_pars_f%from_file )  THEN
     2115!
     2116!--       Horizontal surfaces
     2117          DO  m = 1, surf_lsm_h%ns
     2118             i = surf_lsm_h%i(m)
     2119             j = surf_lsm_h%j(m)
     2120             IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
     2121                surf_lsm_h%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
     2122                surf_lsm_h%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
     2123                surf_lsm_h%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
     2124             ENDIF
     2125          ENDDO
     2126          DO  m = 1, surf_usm_h%ns
     2127             i = surf_usm_h%i(m)
     2128             j = surf_usm_h%j(m)
     2129             IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
     2130                surf_usm_h%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
     2131                surf_usm_h%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
     2132                surf_usm_h%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
     2133             ENDIF
     2134          ENDDO
     2135!
     2136!--       Vertical surfaces
     2137          DO  l = 0, 3
     2138
     2139             ioff = surf_lsm_v(l)%ioff
     2140             joff = surf_lsm_v(l)%joff
     2141             DO  m = 1, surf_lsm_v(l)%ns
     2142                i = surf_lsm_v(l)%i(m) + ioff
     2143                j = surf_lsm_v(l)%j(m) + joff
     2144                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
     2145                   surf_lsm_v(l)%albedo(m,ind_veg_wall) = albedo_pars_f%pars_xy(0,j,i)
     2146                   surf_lsm_v(l)%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
     2147                   surf_lsm_v(l)%albedo(m,ind_wat_win) = albedo_pars_f%pars_xy(0,j,i)
     2148                ENDIF
     2149             ENDDO
     2150
     2151             ioff = surf_usm_v(l)%ioff
     2152             joff = surf_usm_v(l)%joff
     2153             DO  m = 1, surf_usm_v(l)%ns
     2154                i = surf_usm_v(l)%i(m) + ioff
     2155                j = surf_usm_v(l)%j(m) + joff
     2156                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
     2157                   surf_usm_v(l)%albedo(m,ind_veg_wall) = albedo_pars_f%pars_xy(0,j,i)
     2158                   surf_usm_v(l)%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
     2159                   surf_usm_v(l)%albedo(m,ind_wat_win) = albedo_pars_f%pars_xy(0,j,i)
     2160                ENDIF
     2161             ENDDO
     2162          ENDDO
     2163
     2164       ENDIF
     2165!
     2166!--    Read explicit albedo values from building surface pars. If present, they override all less
     2167!--    specific albedo values and force a albedo_type to zero in order to take effect.
     2168       IF ( building_surface_pars_f%from_file )  THEN
     2169          DO  m = 1, surf_usm_h%ns
     2170             i = surf_usm_h%i(m)
     2171             j = surf_usm_h%j(m)
     2172             k = surf_usm_h%k(m)
     2173!
     2174!--          Iterate over surfaces in column, check height and orientation
     2175             DO  is = building_surface_pars_f%index_ji(1,j,i),                                     &
     2176                      building_surface_pars_f%index_ji(2,j,i)
     2177                IF ( building_surface_pars_f%coords(4,is) == -surf_usm_h%koff  .AND.               &
     2178                     building_surface_pars_f%coords(1,is) == k )  THEN
     2179
     2180                   IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=                       &
     2181                        building_surface_pars_f%fill )  THEN
     2182                      surf_usm_h%albedo(m,ind_veg_wall) =                                          &
     2183                               building_surface_pars_f%pars(ind_s_alb_b_wall,is)
     2184                      surf_usm_h%albedo_type(m,ind_veg_wall) = 0
     2185                   ENDIF
     2186
     2187                   IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=                        &
     2188                        building_surface_pars_f%fill )  THEN
     2189                      surf_usm_h%albedo(m,ind_wat_win) =                                           &
     2190                               building_surface_pars_f%pars(ind_s_alb_b_win,is)
     2191                      surf_usm_h%albedo_type(m,ind_wat_win) = 0
     2192                   ENDIF
     2193
     2194                   IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=                      &
     2195                        building_surface_pars_f%fill )  THEN
     2196                      surf_usm_h%albedo(m,ind_pav_green) =                                         &
     2197                               building_surface_pars_f%pars(ind_s_alb_b_green,is)
     2198                      surf_usm_h%albedo_type(m,ind_pav_green) = 0
     2199                   ENDIF
     2200
     2201                   EXIT ! Surface was found and processed
     2202                ENDIF
     2203             ENDDO
     2204          ENDDO
     2205
     2206          DO  l = 0, 3
     2207             DO  m = 1, surf_usm_v(l)%ns
     2208                i = surf_usm_v(l)%i(m)
     2209                j = surf_usm_v(l)%j(m)
     2210                k = surf_usm_v(l)%k(m)
     2211!
     2212!--             Iterate over surfaces in column, check height and orientation
     2213                DO  is = building_surface_pars_f%index_ji(1,j,i),                                  &
     2214                         building_surface_pars_f%index_ji(2,j,i)
     2215                   IF ( building_surface_pars_f%coords(5,is) == -surf_usm_v(l)%joff  .AND.         &
     2216                        building_surface_pars_f%coords(6,is) == -surf_usm_v(l)%ioff  .AND.         &
     2217                        building_surface_pars_f%coords(1,is) == k )  THEN
     2218
     2219                      IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=                    &
     2220                           building_surface_pars_f%fill )  THEN
     2221                         surf_usm_v(l)%albedo(m,ind_veg_wall) =                                    &
     2222                                  building_surface_pars_f%pars(ind_s_alb_b_wall,is)
     2223                         surf_usm_v(l)%albedo_type(m,ind_veg_wall) = 0
     2224                      ENDIF
     2225
     2226                      IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=                     &
     2227                           building_surface_pars_f%fill )  THEN
     2228                         surf_usm_v(l)%albedo(m,ind_wat_win) =                                     &
     2229                                  building_surface_pars_f%pars(ind_s_alb_b_win,is)
     2230                         surf_usm_v(l)%albedo_type(m,ind_wat_win) = 0
     2231                      ENDIF
     2232
     2233                      IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=                   &
     2234                           building_surface_pars_f%fill )  THEN
     2235                         surf_usm_v(l)%albedo(m,ind_pav_green) =                                   &
     2236                                  building_surface_pars_f%pars(ind_s_alb_b_green,is)
     2237                         surf_usm_v(l)%albedo_type(m,ind_pav_green) = 0
     2238                      ENDIF
     2239
     2240                      EXIT ! Surface was found and processed
     2241                   ENDIF
     2242                ENDDO
     2243             ENDDO
     2244          ENDDO
     2245       ENDIF
     2246!
     2247!-- Initialization actions for RRTMG
     2248    ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
     2249#if defined ( __rrtmg )
     2250!
     2251!--    Allocate albedos for short/longwave radiation, horizontal surfaces for wall/green/window
     2252!--    (USM) or vegetation/pavement/water surfaces (LSM).
     2253       ALLOCATE ( surf_lsm_h%aldif(1:surf_lsm_h%ns,0:2)      )
     2254       ALLOCATE ( surf_lsm_h%aldir(1:surf_lsm_h%ns,0:2)      )
     2255       ALLOCATE ( surf_lsm_h%asdif(1:surf_lsm_h%ns,0:2)      )
     2256       ALLOCATE ( surf_lsm_h%asdir(1:surf_lsm_h%ns,0:2)      )
     2257       ALLOCATE ( surf_lsm_h%rrtm_aldif(1:surf_lsm_h%ns,0:2) )
     2258       ALLOCATE ( surf_lsm_h%rrtm_aldir(1:surf_lsm_h%ns,0:2) )
     2259       ALLOCATE ( surf_lsm_h%rrtm_asdif(1:surf_lsm_h%ns,0:2) )
     2260       ALLOCATE ( surf_lsm_h%rrtm_asdir(1:surf_lsm_h%ns,0:2) )
     2261
     2262       ALLOCATE ( surf_usm_h%aldif(1:surf_usm_h%ns,0:2)      )
     2263       ALLOCATE ( surf_usm_h%aldir(1:surf_usm_h%ns,0:2)      )
     2264       ALLOCATE ( surf_usm_h%asdif(1:surf_usm_h%ns,0:2)      )
     2265       ALLOCATE ( surf_usm_h%asdir(1:surf_usm_h%ns,0:2)      )
     2266       ALLOCATE ( surf_usm_h%rrtm_aldif(1:surf_usm_h%ns,0:2) )
     2267       ALLOCATE ( surf_usm_h%rrtm_aldir(1:surf_usm_h%ns,0:2) )
     2268       ALLOCATE ( surf_usm_h%rrtm_asdif(1:surf_usm_h%ns,0:2) )
     2269       ALLOCATE ( surf_usm_h%rrtm_asdir(1:surf_usm_h%ns,0:2) )
     2270
     2271!
     2272!--    Allocate broadband albedo (temporary for the current radiation implementations)
     2273       IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )                                                   &
     2274          ALLOCATE( surf_lsm_h%albedo(1:surf_lsm_h%ns,0:2) )
     2275       IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )                                                   &
     2276          ALLOCATE( surf_usm_h%albedo(1:surf_usm_h%ns,0:2) )
     2277
     2278!
     2279!--    Allocate albedos for short/longwave radiation, vertical surfaces
     2280       DO  l = 0, 3
     2281
     2282          ALLOCATE ( surf_lsm_v(l)%aldif(1:surf_lsm_v(l)%ns,0:2)      )
     2283          ALLOCATE ( surf_lsm_v(l)%aldir(1:surf_lsm_v(l)%ns,0:2)      )
     2284          ALLOCATE ( surf_lsm_v(l)%asdif(1:surf_lsm_v(l)%ns,0:2)      )
     2285          ALLOCATE ( surf_lsm_v(l)%asdir(1:surf_lsm_v(l)%ns,0:2)      )
     2286
     2287          ALLOCATE ( surf_lsm_v(l)%rrtm_aldif(1:surf_lsm_v(l)%ns,0:2) )
     2288          ALLOCATE ( surf_lsm_v(l)%rrtm_aldir(1:surf_lsm_v(l)%ns,0:2) )
     2289          ALLOCATE ( surf_lsm_v(l)%rrtm_asdif(1:surf_lsm_v(l)%ns,0:2) )
     2290          ALLOCATE ( surf_lsm_v(l)%rrtm_asdir(1:surf_lsm_v(l)%ns,0:2) )
     2291
     2292          ALLOCATE ( surf_usm_v(l)%aldif(1:surf_usm_v(l)%ns,0:2)      )
     2293          ALLOCATE ( surf_usm_v(l)%aldir(1:surf_usm_v(l)%ns,0:2)      )
     2294          ALLOCATE ( surf_usm_v(l)%asdif(1:surf_usm_v(l)%ns,0:2)      )
     2295          ALLOCATE ( surf_usm_v(l)%asdir(1:surf_usm_v(l)%ns,0:2)      )
     2296
     2297          ALLOCATE ( surf_usm_v(l)%rrtm_aldif(1:surf_usm_v(l)%ns,0:2) )
     2298          ALLOCATE ( surf_usm_v(l)%rrtm_aldir(1:surf_usm_v(l)%ns,0:2) )
     2299          ALLOCATE ( surf_usm_v(l)%rrtm_asdif(1:surf_usm_v(l)%ns,0:2) )
     2300          ALLOCATE ( surf_usm_v(l)%rrtm_asdir(1:surf_usm_v(l)%ns,0:2) )
     2301!
     2302!--       Allocate broadband albedo (temporary for the current radiation implementations)
     2303          IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )                                           &
     2304             ALLOCATE( surf_lsm_v(l)%albedo(1:surf_lsm_v(l)%ns,0:2) )
     2305          IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )                                           &
     2306             ALLOCATE( surf_usm_v(l)%albedo(1:surf_usm_v(l)%ns,0:2) )
     2307
     2308       ENDDO
     2309!
     2310!--    Level 1 initialization of spectral albedos via namelist paramters. Please note, this case all
     2311!--    surface tiles are initialized the same.
     2312       IF ( surf_lsm_h%ns > 0 )  THEN
     2313          surf_lsm_h%aldif  = albedo_lw_dif
     2314          surf_lsm_h%aldir  = albedo_lw_dir
     2315          surf_lsm_h%asdif  = albedo_sw_dif
     2316          surf_lsm_h%asdir  = albedo_sw_dir
     2317          surf_lsm_h%albedo = albedo_sw_dif
     2318       ENDIF
     2319       IF ( surf_usm_h%ns > 0 )  THEN
     2320          IF ( surf_usm_h%albedo_from_ascii )  THEN
     2321             surf_usm_h%aldif  = surf_usm_h%albedo
     2322             surf_usm_h%aldir  = surf_usm_h%albedo
     2323             surf_usm_h%asdif  = surf_usm_h%albedo
     2324             surf_usm_h%asdir  = surf_usm_h%albedo
     2325          ELSE
     2326             surf_usm_h%aldif  = albedo_lw_dif
     2327             surf_usm_h%aldir  = albedo_lw_dir
     2328             surf_usm_h%asdif  = albedo_sw_dif
     2329             surf_usm_h%asdir  = albedo_sw_dir
     2330             surf_usm_h%albedo = albedo_sw_dif
    17022331          ENDIF
    1703           IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_lw_out_change_0 )  .AND.    &
    1704                      surf_usm_v(l)%ns > 0  )  THEN
    1705              ALLOCATE( surf_usm_v(l)%rad_lw_out_change_0(1:surf_usm_v(l)%ns) )
    1706              surf_usm_v(l)%rad_lw_out_change_0 = 0.0_wp
     2332       ENDIF
     2333
     2334       DO  l = 0, 3
     2335
     2336          IF ( surf_lsm_v(l)%ns > 0 )  THEN
     2337             surf_lsm_v(l)%aldif  = albedo_lw_dif
     2338             surf_lsm_v(l)%aldir  = albedo_lw_dir
     2339             surf_lsm_v(l)%asdif  = albedo_sw_dif
     2340             surf_lsm_v(l)%asdir  = albedo_sw_dir
     2341             surf_lsm_v(l)%albedo = albedo_sw_dif
     2342          ENDIF
     2343
     2344          IF ( surf_usm_v(l)%ns > 0 )  THEN
     2345             IF ( surf_usm_v(l)%albedo_from_ascii )  THEN
     2346                surf_usm_v(l)%aldif  = surf_usm_v(l)%albedo
     2347                surf_usm_v(l)%aldir  = surf_usm_v(l)%albedo
     2348                surf_usm_v(l)%asdif  = surf_usm_v(l)%albedo
     2349                surf_usm_v(l)%asdir  = surf_usm_v(l)%albedo
     2350             ELSE
     2351                surf_usm_v(l)%aldif  = albedo_lw_dif
     2352                surf_usm_v(l)%aldir  = albedo_lw_dir
     2353                surf_usm_v(l)%asdif  = albedo_sw_dif
     2354                surf_usm_v(l)%asdir  = albedo_sw_dir
     2355             ENDIF
    17072356          ENDIF
    17082357       ENDDO
    17092358
    17102359!
    1711 !--    Allocate surface arrays for incoming/outgoing short/longwave radiation
    1712        IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_sw_in )  .AND.                    &
    1713                   surf_lsm_h%ns > 0  )   THEN
    1714           ALLOCATE( surf_lsm_h%rad_sw_in(1:surf_lsm_h%ns)  )
    1715           ALLOCATE( surf_lsm_h%rad_sw_out(1:surf_lsm_h%ns) )
    1716           ALLOCATE( surf_lsm_h%rad_sw_dir(1:surf_lsm_h%ns) )
    1717           ALLOCATE( surf_lsm_h%rad_sw_dif(1:surf_lsm_h%ns) )
    1718           ALLOCATE( surf_lsm_h%rad_sw_ref(1:surf_lsm_h%ns) )
    1719           ALLOCATE( surf_lsm_h%rad_sw_res(1:surf_lsm_h%ns) )
    1720           ALLOCATE( surf_lsm_h%rad_lw_in(1:surf_lsm_h%ns)  )
    1721           ALLOCATE( surf_lsm_h%rad_lw_out(1:surf_lsm_h%ns) )
    1722           ALLOCATE( surf_lsm_h%rad_lw_dif(1:surf_lsm_h%ns) )
    1723           ALLOCATE( surf_lsm_h%rad_lw_ref(1:surf_lsm_h%ns) )
    1724           ALLOCATE( surf_lsm_h%rad_lw_res(1:surf_lsm_h%ns) )
    1725           surf_lsm_h%rad_sw_in  = 0.0_wp
    1726           surf_lsm_h%rad_sw_out = 0.0_wp
    1727           surf_lsm_h%rad_sw_dir = 0.0_wp
    1728           surf_lsm_h%rad_sw_dif = 0.0_wp
    1729           surf_lsm_h%rad_sw_ref = 0.0_wp
    1730           surf_lsm_h%rad_sw_res = 0.0_wp
    1731           surf_lsm_h%rad_lw_in  = 0.0_wp
    1732           surf_lsm_h%rad_lw_out = 0.0_wp
    1733           surf_lsm_h%rad_lw_dif = 0.0_wp
    1734           surf_lsm_h%rad_lw_ref = 0.0_wp
    1735           surf_lsm_h%rad_lw_res = 0.0_wp
    1736        ENDIF
    1737        IF ( .NOT. ALLOCATED ( surf_usm_h%rad_sw_in )  .AND.                    &
    1738                   surf_usm_h%ns > 0  )  THEN
    1739           ALLOCATE( surf_usm_h%rad_sw_in(1:surf_usm_h%ns)  )
    1740           ALLOCATE( surf_usm_h%rad_sw_out(1:surf_usm_h%ns) )
    1741           ALLOCATE( surf_usm_h%rad_sw_dir(1:surf_usm_h%ns) )
    1742           ALLOCATE( surf_usm_h%rad_sw_dif(1:surf_usm_h%ns) )
    1743           ALLOCATE( surf_usm_h%rad_sw_ref(1:surf_usm_h%ns) )
    1744           ALLOCATE( surf_usm_h%rad_sw_res(1:surf_usm_h%ns) )
    1745           ALLOCATE( surf_usm_h%rad_lw_in(1:surf_usm_h%ns)  )
    1746           ALLOCATE( surf_usm_h%rad_lw_out(1:surf_usm_h%ns) )
    1747           ALLOCATE( surf_usm_h%rad_lw_dif(1:surf_usm_h%ns) )
    1748           ALLOCATE( surf_usm_h%rad_lw_ref(1:surf_usm_h%ns) )
    1749           ALLOCATE( surf_usm_h%rad_lw_res(1:surf_usm_h%ns) )
    1750           surf_usm_h%rad_sw_in  = 0.0_wp
    1751           surf_usm_h%rad_sw_out = 0.0_wp
    1752           surf_usm_h%rad_sw_dir = 0.0_wp
    1753           surf_usm_h%rad_sw_dif = 0.0_wp
    1754           surf_usm_h%rad_sw_ref = 0.0_wp
    1755           surf_usm_h%rad_sw_res = 0.0_wp
    1756           surf_usm_h%rad_lw_in  = 0.0_wp
    1757           surf_usm_h%rad_lw_out = 0.0_wp
    1758           surf_usm_h%rad_lw_dif = 0.0_wp
    1759           surf_usm_h%rad_lw_ref = 0.0_wp
    1760           surf_usm_h%rad_lw_res = 0.0_wp
    1761        ENDIF
    1762        DO  l = 0, 3
    1763           IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_sw_in )  .AND.              &
    1764                      surf_lsm_v(l)%ns > 0  )  THEN
    1765              ALLOCATE( surf_lsm_v(l)%rad_sw_in(1:surf_lsm_v(l)%ns)  )
    1766              ALLOCATE( surf_lsm_v(l)%rad_sw_out(1:surf_lsm_v(l)%ns) )
    1767              ALLOCATE( surf_lsm_v(l)%rad_sw_dir(1:surf_lsm_v(l)%ns) )
    1768              ALLOCATE( surf_lsm_v(l)%rad_sw_dif(1:surf_lsm_v(l)%ns) )
    1769              ALLOCATE( surf_lsm_v(l)%rad_sw_ref(1:surf_lsm_v(l)%ns) )
    1770              ALLOCATE( surf_lsm_v(l)%rad_sw_res(1:surf_lsm_v(l)%ns) )
    1771 
    1772              ALLOCATE( surf_lsm_v(l)%rad_lw_in(1:surf_lsm_v(l)%ns)  )
    1773              ALLOCATE( surf_lsm_v(l)%rad_lw_out(1:surf_lsm_v(l)%ns) )
    1774              ALLOCATE( surf_lsm_v(l)%rad_lw_dif(1:surf_lsm_v(l)%ns) )
    1775              ALLOCATE( surf_lsm_v(l)%rad_lw_ref(1:surf_lsm_v(l)%ns) )
    1776              ALLOCATE( surf_lsm_v(l)%rad_lw_res(1:surf_lsm_v(l)%ns) )
    1777 
    1778              surf_lsm_v(l)%rad_sw_in  = 0.0_wp
    1779              surf_lsm_v(l)%rad_sw_out = 0.0_wp
    1780              surf_lsm_v(l)%rad_sw_dir = 0.0_wp
    1781              surf_lsm_v(l)%rad_sw_dif = 0.0_wp
    1782              surf_lsm_v(l)%rad_sw_ref = 0.0_wp
    1783              surf_lsm_v(l)%rad_sw_res = 0.0_wp
    1784 
    1785              surf_lsm_v(l)%rad_lw_in  = 0.0_wp
    1786              surf_lsm_v(l)%rad_lw_out = 0.0_wp
    1787              surf_lsm_v(l)%rad_lw_dif = 0.0_wp
    1788              surf_lsm_v(l)%rad_lw_ref = 0.0_wp
    1789              surf_lsm_v(l)%rad_lw_res = 0.0_wp
    1790           ENDIF
    1791           IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_sw_in )  .AND.              &
    1792                      surf_usm_v(l)%ns > 0  )  THEN
    1793              ALLOCATE( surf_usm_v(l)%rad_sw_in(1:surf_usm_v(l)%ns)  )
    1794              ALLOCATE( surf_usm_v(l)%rad_sw_out(1:surf_usm_v(l)%ns) )
    1795              ALLOCATE( surf_usm_v(l)%rad_sw_dir(1:surf_usm_v(l)%ns) )
    1796              ALLOCATE( surf_usm_v(l)%rad_sw_dif(1:surf_usm_v(l)%ns) )
    1797              ALLOCATE( surf_usm_v(l)%rad_sw_ref(1:surf_usm_v(l)%ns) )
    1798              ALLOCATE( surf_usm_v(l)%rad_sw_res(1:surf_usm_v(l)%ns) )
    1799              ALLOCATE( surf_usm_v(l)%rad_lw_in(1:surf_usm_v(l)%ns)  )
    1800              ALLOCATE( surf_usm_v(l)%rad_lw_out(1:surf_usm_v(l)%ns) )
    1801              ALLOCATE( surf_usm_v(l)%rad_lw_dif(1:surf_usm_v(l)%ns) )
    1802              ALLOCATE( surf_usm_v(l)%rad_lw_ref(1:surf_usm_v(l)%ns) )
    1803              ALLOCATE( surf_usm_v(l)%rad_lw_res(1:surf_usm_v(l)%ns) )
    1804              surf_usm_v(l)%rad_sw_in  = 0.0_wp
    1805              surf_usm_v(l)%rad_sw_out = 0.0_wp
    1806              surf_usm_v(l)%rad_sw_dir = 0.0_wp
    1807              surf_usm_v(l)%rad_sw_dif = 0.0_wp
    1808              surf_usm_v(l)%rad_sw_ref = 0.0_wp
    1809              surf_usm_v(l)%rad_sw_res = 0.0_wp
    1810              surf_usm_v(l)%rad_lw_in  = 0.0_wp
    1811              surf_usm_v(l)%rad_lw_out = 0.0_wp
    1812              surf_usm_v(l)%rad_lw_dif = 0.0_wp
    1813              surf_usm_v(l)%rad_lw_ref = 0.0_wp
    1814              surf_usm_v(l)%rad_lw_res = 0.0_wp
     2360!--    Level 2 initialization of spectral albedos via albedo_type.
     2361!--    Please note, for natural- and urban-type surfaces, a tile approach is applied so that the
     2362!--    resulting albedo is calculated via the weighted average of respective surface fractions.
     2363       DO  m = 1, surf_lsm_h%ns
     2364!
     2365!--       Spectral albedos for vegetation/pavement/water surfaces
     2366          DO  ind_type = 0, 2
     2367             IF ( surf_lsm_h%albedo_type(m,ind_type) /= 0 )  THEN
     2368                surf_lsm_h%aldif(m,ind_type) = albedo_pars(1,surf_lsm_h%albedo_type(m,ind_type))
     2369                surf_lsm_h%asdif(m,ind_type) = albedo_pars(2,surf_lsm_h%albedo_type(m,ind_type))
     2370                surf_lsm_h%aldir(m,ind_type) = albedo_pars(1,surf_lsm_h%albedo_type(m,ind_type))
     2371                surf_lsm_h%asdir(m,ind_type) = albedo_pars(2,surf_lsm_h%albedo_type(m,ind_type))
     2372                surf_lsm_h%albedo(m,ind_type) = albedo_pars(0,surf_lsm_h%albedo_type(m,ind_type))
     2373             ENDIF
     2374          ENDDO
     2375
     2376       ENDDO
     2377!
     2378!--    For urban surface only if albedo has not already been initialized in the urban-surface model
     2379!--    via the ASCII file.
     2380       IF ( .NOT. surf_usm_h%albedo_from_ascii )  THEN
     2381          DO  m = 1, surf_usm_h%ns
     2382!
     2383!--          Spectral albedos for wall/green/window surfaces
     2384             DO  ind_type = 0, 2
     2385                IF ( surf_usm_h%albedo_type(m,ind_type) /= 0 )  THEN
     2386                   surf_usm_h%aldif(m,ind_type) = albedo_pars(1,surf_usm_h%albedo_type(m,ind_type))
     2387                   surf_usm_h%asdif(m,ind_type) = albedo_pars(2,surf_usm_h%albedo_type(m,ind_type))
     2388                   surf_usm_h%aldir(m,ind_type) = albedo_pars(1,surf_usm_h%albedo_type(m,ind_type))
     2389                   surf_usm_h%asdir(m,ind_type) = albedo_pars(2,surf_usm_h%albedo_type(m,ind_type))
     2390                   surf_usm_h%albedo(m,ind_type) = albedo_pars(0,surf_usm_h%albedo_type(m,ind_type))
     2391                ENDIF
     2392             ENDDO
     2393
     2394          ENDDO
     2395       ENDIF
     2396
     2397       DO l = 0, 3
     2398
     2399          DO  m = 1, surf_lsm_v(l)%ns
     2400!
     2401!--          Spectral albedos for vegetation/pavement/water surfaces
     2402             DO  ind_type = 0, 2
     2403                IF ( surf_lsm_v(l)%albedo_type(m,ind_type) /= 0 )  THEN
     2404                   surf_lsm_v(l)%aldif(m,ind_type) =                                               &
     2405                         albedo_pars(1,surf_lsm_v(l)%albedo_type(m,ind_type))
     2406                   surf_lsm_v(l)%asdif(m,ind_type) =                                               &
     2407                         albedo_pars(2,surf_lsm_v(l)%albedo_type(m,ind_type))
     2408                   surf_lsm_v(l)%aldir(m,ind_type) =                                               &
     2409                         albedo_pars(1,surf_lsm_v(l)%albedo_type(m,ind_type))
     2410                   surf_lsm_v(l)%asdir(m,ind_type) =                                               &
     2411                         albedo_pars(2,surf_lsm_v(l)%albedo_type(m,ind_type))
     2412                   surf_lsm_v(l)%albedo(m,ind_type) =                                              &
     2413                         albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_type))
     2414                ENDIF
     2415             ENDDO
     2416          ENDDO
     2417!
     2418!--       For urban surface only if albedo has not already been initialized in the urban-surface
     2419!--       model via the ASCII file.
     2420          IF ( .NOT. surf_usm_v(l)%albedo_from_ascii )  THEN
     2421             DO  m = 1, surf_usm_v(l)%ns
     2422!
     2423!--             Spectral albedos for wall/green/window surfaces
     2424                DO  ind_type = 0, 2
     2425                   IF ( surf_usm_v(l)%albedo_type(m,ind_type) /= 0 )  THEN
     2426                      surf_usm_v(l)%aldif(m,ind_type) =                                            &
     2427                         albedo_pars(1,surf_usm_v(l)%albedo_type(m,ind_type))
     2428                      surf_usm_v(l)%asdif(m,ind_type) =                                            &
     2429                         albedo_pars(2,surf_usm_v(l)%albedo_type(m,ind_type))
     2430                      surf_usm_v(l)%aldir(m,ind_type) =                                            &
     2431                         albedo_pars(1,surf_usm_v(l)%albedo_type(m,ind_type))
     2432                      surf_usm_v(l)%asdir(m,ind_type) =                                            &
     2433                         albedo_pars(2,surf_usm_v(l)%albedo_type(m,ind_type))
     2434                      surf_usm_v(l)%albedo(m,ind_type) =                                           &
     2435                         albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_type))
     2436                   ENDIF
     2437                ENDDO
     2438
     2439             ENDDO
    18152440          ENDIF
    18162441       ENDDO
    18172442!
    1818 !--    Fix net radiation in case of radiation_scheme = 'constant'
    1819        IF ( radiation_scheme == 'constant' )  THEN
    1820           IF ( ALLOCATED( surf_lsm_h%rad_net ) )                               &
    1821              surf_lsm_h%rad_net    = net_radiation
    1822           IF ( ALLOCATED( surf_usm_h%rad_net ) )                               &
    1823              surf_usm_h%rad_net    = net_radiation
    1824 !
    1825 !--       Todo: weight with inclination angle
    1826           DO  l = 0, 3
    1827              IF ( ALLOCATED( surf_lsm_v(l)%rad_net ) )                         &
    1828                 surf_lsm_v(l)%rad_net = net_radiation
    1829              IF ( ALLOCATED( surf_usm_v(l)%rad_net ) )                         &
    1830                 surf_usm_v(l)%rad_net = net_radiation
    1831           ENDDO
    1832 !          radiation = .FALSE.
    1833 !
    1834 !--    Calculate orbital constants
    1835        ELSE
    1836           decl_1 = SIN(23.45_wp * pi / 180.0_wp)
    1837           decl_2 = 2.0_wp * pi / 365.0_wp
    1838           decl_3 = decl_2 * 81.0_wp
    1839           lat    = latitude * pi / 180.0_wp
    1840           lon    = longitude * pi / 180.0_wp
    1841        ENDIF
    1842 
    1843        IF ( radiation_scheme == 'clear-sky'  .OR.                              &
    1844             radiation_scheme == 'constant'   .OR.                              &
    1845             radiation_scheme == 'external' )  THEN
    1846 !
    1847 !--       Allocate arrays for incoming/outgoing short/longwave radiation
    1848           IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
    1849              ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
    1850              rad_sw_in = 0.0_wp
    1851           ENDIF
    1852           IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
    1853              ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
    1854              rad_sw_out = 0.0_wp
    1855           ENDIF
    1856 
    1857           IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
    1858              ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
    1859              rad_lw_in = 0.0_wp
    1860           ENDIF
    1861           IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
    1862              ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
    1863              rad_lw_out = 0.0_wp
    1864           ENDIF
    1865 
    1866 !
    1867 !--       Allocate average arrays for incoming/outgoing short/longwave radiation
    1868           IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
    1869              ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
    1870           ENDIF
    1871           IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
    1872              ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
    1873           ENDIF
    1874 
    1875           IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
    1876              ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
    1877           ENDIF
    1878           IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
    1879              ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
    1880           ENDIF
    1881 !
    1882 !--       Allocate arrays for broadband albedo, and level 1 initialization
    1883 !--       via namelist paramter, unless not already allocated.
    1884           IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )  THEN
    1885              ALLOCATE( surf_lsm_h%albedo(1:surf_lsm_h%ns,0:2)     )
    1886              surf_lsm_h%albedo    = albedo
    1887           ENDIF
    1888           IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )  THEN
    1889              ALLOCATE( surf_usm_h%albedo(1:surf_usm_h%ns,0:2)     )
    1890              surf_usm_h%albedo    = albedo
    1891           ENDIF
    1892 
    1893           DO  l = 0, 3
    1894              IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )  THEN
    1895                 ALLOCATE( surf_lsm_v(l)%albedo(1:surf_lsm_v(l)%ns,0:2) )
    1896                 surf_lsm_v(l)%albedo = albedo
    1897              ENDIF
    1898              IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )  THEN
    1899                 ALLOCATE( surf_usm_v(l)%albedo(1:surf_usm_v(l)%ns,0:2) )
    1900                 surf_usm_v(l)%albedo = albedo
    1901              ENDIF
    1902           ENDDO
    1903 !
    1904 !--       Level 2 initialization of broadband albedo via given albedo_type.
    1905 !--       Only if albedo_type is non-zero. In case of urban surface and
    1906 !--       input data is read from ASCII file, albedo_type will be zero, so that
    1907 !--       albedo won't be overwritten.
     2443!--    Level 3 initialization at grid points where albedo type is zero.
     2444!--    In this case, spectral albedos are taken from file if available
     2445       IF ( albedo_pars_f%from_file )  THEN
     2446!
     2447!--       Horizontal
    19082448          DO  m = 1, surf_lsm_h%ns
    1909              IF ( surf_lsm_h%albedo_type(m,ind_veg_wall) /= 0 )                &
    1910                 surf_lsm_h%albedo(m,ind_veg_wall) =                            &
    1911                            albedo_pars(0,surf_lsm_h%albedo_type(m,ind_veg_wall))
    1912              IF ( surf_lsm_h%albedo_type(m,ind_pav_green) /= 0 )               &
    1913                 surf_lsm_h%albedo(m,ind_pav_green) =                           &
    1914                            albedo_pars(0,surf_lsm_h%albedo_type(m,ind_pav_green))
    1915              IF ( surf_lsm_h%albedo_type(m,ind_wat_win) /= 0 )                 &
    1916                 surf_lsm_h%albedo(m,ind_wat_win) =                             &
    1917                            albedo_pars(0,surf_lsm_h%albedo_type(m,ind_wat_win))
    1918           ENDDO
    1919           DO  m = 1, surf_usm_h%ns
    1920              IF ( surf_usm_h%albedo_type(m,ind_veg_wall) /= 0 )                &
    1921                 surf_usm_h%albedo(m,ind_veg_wall) =                            &
    1922                            albedo_pars(0,surf_usm_h%albedo_type(m,ind_veg_wall))
    1923              IF ( surf_usm_h%albedo_type(m,ind_pav_green) /= 0 )               &
    1924                 surf_usm_h%albedo(m,ind_pav_green) =                           &
    1925                            albedo_pars(0,surf_usm_h%albedo_type(m,ind_pav_green))
    1926              IF ( surf_usm_h%albedo_type(m,ind_wat_win) /= 0 )                 &
    1927                 surf_usm_h%albedo(m,ind_wat_win) =                             &
    1928                            albedo_pars(0,surf_usm_h%albedo_type(m,ind_wat_win))
    1929           ENDDO
    1930 
    1931           DO  l = 0, 3
    1932              DO  m = 1, surf_lsm_v(l)%ns
    1933                 IF ( surf_lsm_v(l)%albedo_type(m,ind_veg_wall) /= 0 )          &
    1934                    surf_lsm_v(l)%albedo(m,ind_veg_wall) =                      &
    1935                         albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_veg_wall))
    1936                 IF ( surf_lsm_v(l)%albedo_type(m,ind_pav_green) /= 0 )         &
    1937                    surf_lsm_v(l)%albedo(m,ind_pav_green) =                     &
    1938                         albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_pav_green))
    1939                 IF ( surf_lsm_v(l)%albedo_type(m,ind_wat_win) /= 0 )           &
    1940                    surf_lsm_v(l)%albedo(m,ind_wat_win) =                       &
    1941                         albedo_pars(0,surf_lsm_v(l)%albedo_type(m,ind_wat_win))
    1942              ENDDO
    1943              DO  m = 1, surf_usm_v(l)%ns
    1944                 IF ( surf_usm_v(l)%albedo_type(m,ind_veg_wall) /= 0 )          &
    1945                    surf_usm_v(l)%albedo(m,ind_veg_wall) =                      &
    1946                         albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_veg_wall))
    1947                 IF ( surf_usm_v(l)%albedo_type(m,ind_pav_green) /= 0 )         &
    1948                    surf_usm_v(l)%albedo(m,ind_pav_green) =                     &
    1949                         albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_pav_green))
    1950                 IF ( surf_usm_v(l)%albedo_type(m,ind_wat_win) /= 0 )           &
    1951                    surf_usm_v(l)%albedo(m,ind_wat_win) =                       &
    1952                         albedo_pars(0,surf_usm_v(l)%albedo_type(m,ind_wat_win))
     2449             i = surf_lsm_h%i(m)
     2450             j = surf_lsm_h%j(m)
     2451!
     2452!--          Spectral albedos for vegetation/pavement/water surfaces
     2453             DO  ind_type = 0, 2
     2454                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )                          &
     2455                   surf_lsm_h%albedo(m,ind_type) =  albedo_pars_f%pars_xy(0,j,i)
     2456                IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )                          &
     2457                   surf_lsm_h%aldir(m,ind_type) = albedo_pars_f%pars_xy(1,j,i)
     2458                IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )                          &
     2459                   surf_lsm_h%aldif(m,ind_type) = albedo_pars_f%pars_xy(1,j,i)
     2460                IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )                          &
     2461                   surf_lsm_h%asdir(m,ind_type) = albedo_pars_f%pars_xy(2,j,i)
     2462                IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )                          &
     2463                   surf_lsm_h%asdif(m,ind_type) = albedo_pars_f%pars_xy(2,j,i)
    19532464             ENDDO
    19542465          ENDDO
    1955 
    1956 !
    1957 !--       Level 3 initialization at grid points where albedo type is zero.
    1958 !--       This case, albedo is taken from file. In case of constant radiation
    1959 !--       or clear sky, only broadband albedo is given.
    1960           IF ( albedo_pars_f%from_file )  THEN
    1961 !
    1962 !--          Horizontal surfaces
    1963              DO  m = 1, surf_lsm_h%ns
    1964                 i = surf_lsm_h%i(m)
    1965                 j = surf_lsm_h%j(m)
    1966                 IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
    1967                    surf_lsm_h%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
    1968                    surf_lsm_h%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
    1969                    surf_lsm_h%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
    1970                 ENDIF
    1971              ENDDO
     2466!
     2467!--       For urban surface only if albedo has not been already initialized in the urban-surface
     2468!--       model via the ASCII file.
     2469          IF ( .NOT. surf_usm_h%albedo_from_ascii )  THEN
    19722470             DO  m = 1, surf_usm_h%ns
    19732471                i = surf_usm_h%i(m)
    19742472                j = surf_usm_h%j(m)
    1975                 IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
    1976                    surf_usm_h%albedo(m,ind_veg_wall)  = albedo_pars_f%pars_xy(0,j,i)
    1977                    surf_usm_h%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
    1978                    surf_usm_h%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
     2473!
     2474!--             Broadband albedos for wall/green/window surfaces
     2475                DO  ind_type = 0, 2
     2476                   IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )                       &
     2477                      surf_usm_h%albedo(m,ind_type) = albedo_pars_f%pars_xy(0,j,i)
     2478                ENDDO
     2479!
     2480!--             Spectral albedos especially for building wall surfaces
     2481                IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )  THEN
     2482                   surf_usm_h%aldir(m,ind_veg_wall) = albedo_pars_f%pars_xy(1,j,i)
     2483                   surf_usm_h%aldif(m,ind_veg_wall) = albedo_pars_f%pars_xy(1,j,i)
     2484                ENDIF
     2485                IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )  THEN
     2486                   surf_usm_h%asdir(m,ind_veg_wall) = albedo_pars_f%pars_xy(2,j,i)
     2487                   surf_usm_h%asdif(m,ind_veg_wall) = albedo_pars_f%pars_xy(2,j,i)
     2488                ENDIF
     2489!
     2490!--             Spectral albedos especially for building green surfaces
     2491                IF ( albedo_pars_f%pars_xy(3,j,i) /= albedo_pars_f%fill )  THEN
     2492                   surf_usm_h%aldir(m,ind_pav_green) = albedo_pars_f%pars_xy(3,j,i)
     2493                   surf_usm_h%aldif(m,ind_pav_green) = albedo_pars_f%pars_xy(3,j,i)
     2494                ENDIF
     2495                IF ( albedo_pars_f%pars_xy(4,j,i) /= albedo_pars_f%fill )  THEN
     2496                   surf_usm_h%asdir(m,ind_pav_green) = albedo_pars_f%pars_xy(4,j,i)
     2497                   surf_usm_h%asdif(m,ind_pav_green) = albedo_pars_f%pars_xy(4,j,i)
     2498                ENDIF
     2499!
     2500!--             Spectral albedos especially for building window surfaces
     2501                IF ( albedo_pars_f%pars_xy(5,j,i) /= albedo_pars_f%fill )  THEN
     2502                   surf_usm_h%aldir(m,ind_wat_win) = albedo_pars_f%pars_xy(5,j,i)
     2503                   surf_usm_h%aldif(m,ind_wat_win) = albedo_pars_f%pars_xy(5,j,i)
     2504                ENDIF
     2505                IF ( albedo_pars_f%pars_xy(6,j,i) /= albedo_pars_f%fill )  THEN
     2506                   surf_usm_h%asdir(m,ind_wat_win) = albedo_pars_f%pars_xy(6,j,i)
     2507                   surf_usm_h%asdif(m,ind_wat_win) = albedo_pars_f%pars_xy(6,j,i)
     2508                ENDIF
     2509
     2510             ENDDO
     2511          ENDIF
     2512!
     2513!--       Vertical
     2514          DO  l = 0, 3
     2515             ioff = surf_lsm_v(l)%ioff
     2516             joff = surf_lsm_v(l)%joff
     2517
     2518             DO  m = 1, surf_lsm_v(l)%ns
     2519                i = surf_lsm_v(l)%i(m)
     2520                j = surf_lsm_v(l)%j(m)
     2521!
     2522!--             Spectral albedos for vegetation/pavement/water surfaces
     2523                DO  ind_type = 0, 2
     2524                   IF ( albedo_pars_f%pars_xy(0,j+joff,i+ioff) /= albedo_pars_f%fill )             &
     2525                      surf_lsm_v(l)%albedo(m,ind_type) = albedo_pars_f%pars_xy(0,j+joff,i+ioff)
     2526                   IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /= albedo_pars_f%fill )             &
     2527                      surf_lsm_v(l)%aldir(m,ind_type) = albedo_pars_f%pars_xy(1,j+joff,i+ioff)
     2528                   IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /= albedo_pars_f%fill )             &
     2529                      surf_lsm_v(l)%aldif(m,ind_type) = albedo_pars_f%pars_xy(1,j+joff,i+ioff)
     2530                   IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /= albedo_pars_f%fill )             &
     2531                      surf_lsm_v(l)%asdir(m,ind_type) = albedo_pars_f%pars_xy(2,j+joff,i+ioff)
     2532                   IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /= albedo_pars_f%fill )             &
     2533                      surf_lsm_v(l)%asdif(m,ind_type) = albedo_pars_f%pars_xy(2,j+joff,i+ioff)
     2534                ENDDO
     2535             ENDDO
     2536!
     2537!--          For urban surface only if albedo has not already been initialized in the urban-surface
     2538!--          model via the ASCII file.
     2539             IF ( .NOT. surf_usm_v(l)%albedo_from_ascii )  THEN
     2540                ioff = surf_usm_v(l)%ioff
     2541                joff = surf_usm_v(l)%joff
     2542
     2543                DO  m = 1, surf_usm_v(l)%ns
     2544                   i = surf_usm_v(l)%i(m)
     2545                   j = surf_usm_v(l)%j(m)
     2546!
     2547!--                Broadband albedos for wall/green/window surfaces
     2548                   DO  ind_type = 0, 2
     2549                      IF ( albedo_pars_f%pars_xy(0,j+joff,i+ioff) /= albedo_pars_f%fill )          &
     2550                         surf_usm_v(l)%albedo(m,ind_type) = albedo_pars_f%pars_xy(0,j+joff,i+ioff)
     2551                   ENDDO
     2552!
     2553!--                Spectral albedos especially for building wall surfaces
     2554                   IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /= albedo_pars_f%fill )  THEN
     2555                      surf_usm_v(l)%aldir(m,ind_veg_wall) = albedo_pars_f%pars_xy(1,j+joff,i+ioff)
     2556                      surf_usm_v(l)%aldif(m,ind_veg_wall) = albedo_pars_f%pars_xy(1,j+joff,i+ioff)
     2557                   ENDIF
     2558                   IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /= albedo_pars_f%fill )  THEN
     2559                      surf_usm_v(l)%asdir(m,ind_veg_wall) = albedo_pars_f%pars_xy(2,j+joff,i+ioff)
     2560                      surf_usm_v(l)%asdif(m,ind_veg_wall) = albedo_pars_f%pars_xy(2,j+joff,i+ioff)
     2561                   ENDIF
     2562!
     2563!--                Spectral albedos especially for building green surfaces
     2564                   IF ( albedo_pars_f%pars_xy(3,j+joff,i+ioff) /= albedo_pars_f%fill )  THEN
     2565                      surf_usm_v(l)%aldir(m,ind_pav_green) = albedo_pars_f%pars_xy(3,j+joff,i+ioff)
     2566                      surf_usm_v(l)%aldif(m,ind_pav_green) = albedo_pars_f%pars_xy(3,j+joff,i+ioff)
     2567                   ENDIF
     2568                   IF ( albedo_pars_f%pars_xy(4,j+joff,i+ioff) /= albedo_pars_f%fill )  THEN
     2569                      surf_usm_v(l)%asdir(m,ind_pav_green) = albedo_pars_f%pars_xy(4,j+joff,i+ioff)
     2570                      surf_usm_v(l)%asdif(m,ind_pav_green) = albedo_pars_f%pars_xy(4,j+joff,i+ioff)
     2571                   ENDIF
     2572!
     2573!--                Spectral albedos especially for building window surfaces
     2574                   IF ( albedo_pars_f%pars_xy(5,j+joff,i+ioff) /= albedo_pars_f%fill )  THEN
     2575                      surf_usm_v(l)%aldir(m,ind_wat_win) = albedo_pars_f%pars_xy(5,j+joff,i+ioff)
     2576                      surf_usm_v(l)%aldif(m,ind_wat_win) = albedo_pars_f%pars_xy(5,j+joff,i+ioff)
     2577                   ENDIF
     2578                   IF ( albedo_pars_f%pars_xy(6,j+joff,i+ioff) /= albedo_pars_f%fill )  THEN
     2579                      surf_usm_v(l)%asdir(m,ind_wat_win) = albedo_pars_f%pars_xy(6,j+joff,i+ioff)
     2580                      surf_usm_v(l)%asdif(m,ind_wat_win) = albedo_pars_f%pars_xy(6,j+joff,i+ioff)
     2581                   ENDIF
     2582                ENDDO
     2583             ENDIF
     2584          ENDDO
     2585
     2586       ENDIF
     2587!
     2588!--    Read explicit albedo values from building surface pars. If present they override all less
     2589!--    specific albedo values and force an albedo_type to zero in order to take effect.
     2590       IF ( building_surface_pars_f%from_file )  THEN
     2591          DO  m = 1, surf_usm_h%ns
     2592             i = surf_usm_h%i(m)
     2593             j = surf_usm_h%j(m)
     2594             k = surf_usm_h%k(m)
     2595!
     2596!--          Iterate over surfaces in column, check height and orientation
     2597             DO  is = building_surface_pars_f%index_ji(1,j,i),                                     &
     2598                      building_surface_pars_f%index_ji(2,j,i)
     2599                IF ( building_surface_pars_f%coords(4,is) == -surf_usm_h%koff  .AND.               &
     2600                     building_surface_pars_f%coords(1,is) == k )  THEN
     2601
     2602                   IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=                       &
     2603                        building_surface_pars_f%fill )  THEN
     2604                      surf_usm_h%albedo(m,ind_veg_wall) =                                          &
     2605                               building_surface_pars_f%pars(ind_s_alb_b_wall,is)
     2606                      surf_usm_h%albedo_type(m,ind_veg_wall) = 0
     2607                   ENDIF
     2608
     2609                   IF ( building_surface_pars_f%pars(ind_s_alb_l_wall,is) /=                       &
     2610                        building_surface_pars_f%fill )  THEN
     2611                      surf_usm_h%aldir(m,ind_veg_wall) =                                           &
     2612                               building_surface_pars_f%pars(ind_s_alb_l_wall,is)
     2613                      surf_usm_h%aldif(m,ind_veg_wall) =                                           &
     2614                               building_surface_pars_f%pars(ind_s_alb_l_wall,is)
     2615                      surf_usm_h%albedo_type(m,ind_veg_wall) = 0
     2616                   ENDIF
     2617
     2618                   IF ( building_surface_pars_f%pars(ind_s_alb_s_wall,is) /=                       &
     2619                        building_surface_pars_f%fill )  THEN
     2620                      surf_usm_h%asdir(m,ind_veg_wall) =                                           &
     2621                               building_surface_pars_f%pars(ind_s_alb_s_wall,is)
     2622                      surf_usm_h%asdif(m,ind_veg_wall) =                                           &
     2623                               building_surface_pars_f%pars(ind_s_alb_s_wall,is)
     2624                      surf_usm_h%albedo_type(m,ind_veg_wall) = 0
     2625                   ENDIF
     2626
     2627                   IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=                        &
     2628                        building_surface_pars_f%fill )  THEN
     2629                      surf_usm_h%albedo(m,ind_wat_win) =                                           &
     2630                               building_surface_pars_f%pars(ind_s_alb_b_win,is)
     2631                      surf_usm_h%albedo_type(m,ind_wat_win) = 0
     2632                   ENDIF
     2633
     2634                   IF ( building_surface_pars_f%pars(ind_s_alb_l_win,is) /=                        &
     2635                        building_surface_pars_f%fill )  THEN
     2636                      surf_usm_h%aldir(m,ind_wat_win) =                                            &
     2637                               building_surface_pars_f%pars(ind_s_alb_l_win,is)
     2638                      surf_usm_h%aldif(m,ind_wat_win) =                                            &
     2639                               building_surface_pars_f%pars(ind_s_alb_l_win,is)
     2640                      surf_usm_h%albedo_type(m,ind_wat_win) = 0
     2641                   ENDIF
     2642
     2643                   IF ( building_surface_pars_f%pars(ind_s_alb_s_win,is) /=                        &
     2644                        building_surface_pars_f%fill )  THEN
     2645                      surf_usm_h%asdir(m,ind_wat_win) =                                            &
     2646                               building_surface_pars_f%pars(ind_s_alb_s_win,is)
     2647                      surf_usm_h%asdif(m,ind_wat_win) =                                            &
     2648                               building_surface_pars_f%pars(ind_s_alb_s_win,is)
     2649                      surf_usm_h%albedo_type(m,ind_wat_win) = 0
     2650                   ENDIF
     2651
     2652                   IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=                      &
     2653                        building_surface_pars_f%fill )  THEN
     2654                      surf_usm_h%albedo(m,ind_pav_green) =                                         &
     2655                               building_surface_pars_f%pars(ind_s_alb_b_green,is)
     2656                      surf_usm_h%albedo_type(m,ind_pav_green) = 0
     2657                   ENDIF
     2658
     2659                   IF ( building_surface_pars_f%pars(ind_s_alb_l_green,is) /=                      &
     2660                        building_surface_pars_f%fill )  THEN
     2661                      surf_usm_h%aldir(m,ind_pav_green) =                                          &
     2662                               building_surface_pars_f%pars(ind_s_alb_l_green,is)
     2663                      surf_usm_h%aldif(m,ind_pav_green) =                                          &
     2664                               building_surface_pars_f%pars(ind_s_alb_l_green,is)
     2665                      surf_usm_h%albedo_type(m,ind_pav_green) = 0
     2666                   ENDIF
     2667
     2668                   IF ( building_surface_pars_f%pars(ind_s_alb_s_green,is) /=                      &
     2669                        building_surface_pars_f%fill )  THEN
     2670                      surf_usm_h%asdir(m,ind_pav_green) =                                          &
     2671                               building_surface_pars_f%pars(ind_s_alb_s_green,is)
     2672                      surf_usm_h%asdif(m,ind_pav_green) =                                          &
     2673                               building_surface_pars_f%pars(ind_s_alb_s_green,is)
     2674                      surf_usm_h%albedo_type(m,ind_pav_green) = 0
     2675                   ENDIF
     2676
     2677                   EXIT ! Surface was found and processed
    19792678                ENDIF
    19802679             ENDDO
    1981 !
    1982 !--          Vertical surfaces
    1983              DO  l = 0, 3
    1984 
    1985                 ioff = surf_lsm_v(l)%ioff
    1986                 joff = surf_lsm_v(l)%joff
    1987                 DO  m = 1, surf_lsm_v(l)%ns
    1988                    i = surf_lsm_v(l)%i(m) + ioff
    1989                    j = surf_lsm_v(l)%j(m) + joff
    1990                    IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
    1991                       surf_lsm_v(l)%albedo(m,ind_veg_wall) = albedo_pars_f%pars_xy(0,j,i)
    1992                       surf_lsm_v(l)%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
    1993                       surf_lsm_v(l)%albedo(m,ind_wat_win) = albedo_pars_f%pars_xy(0,j,i)
    1994                    ENDIF
    1995                 ENDDO
    1996 
    1997                 ioff = surf_usm_v(l)%ioff
    1998                 joff = surf_usm_v(l)%joff
    1999                 DO  m = 1, surf_usm_v(l)%ns
    2000                    i = surf_usm_v(l)%i(m) + ioff
    2001                    j = surf_usm_v(l)%j(m) + joff
    2002                    IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
    2003                       surf_usm_v(l)%albedo(m,ind_veg_wall) = albedo_pars_f%pars_xy(0,j,i)
    2004                       surf_usm_v(l)%albedo(m,ind_pav_green) = albedo_pars_f%pars_xy(0,j,i)
    2005                       surf_usm_v(l)%albedo(m,ind_wat_win) = albedo_pars_f%pars_xy(0,j,i)
     2680          ENDDO
     2681
     2682          DO  l = 0, 3
     2683             DO  m = 1, surf_usm_v(l)%ns
     2684                i = surf_usm_v(l)%i(m)
     2685                j = surf_usm_v(l)%j(m)
     2686                k = surf_usm_v(l)%k(m)
     2687!
     2688!--             Iterate over surfaces in column, check height and orientation
     2689                DO  is = building_surface_pars_f%index_ji(1,j,i),                                  &
     2690                         building_surface_pars_f%index_ji(2,j,i)
     2691                   IF ( building_surface_pars_f%coords(5,is) == -surf_usm_v(l)%joff  .AND.         &
     2692                        building_surface_pars_f%coords(6,is) == -surf_usm_v(l)%ioff  .AND.         &
     2693                        building_surface_pars_f%coords(1,is) == k )  THEN
     2694
     2695                      IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=                    &
     2696                           building_surface_pars_f%fill )  THEN
     2697                         surf_usm_v(l)%albedo(m,ind_veg_wall) =                                    &
     2698                                  building_surface_pars_f%pars(ind_s_alb_b_wall,is)
     2699                         surf_usm_v(l)%albedo_type(m,ind_veg_wall) = 0
     2700                      ENDIF
     2701
     2702                      IF ( building_surface_pars_f%pars(ind_s_alb_l_wall,is) /=                    &
     2703                           building_surface_pars_f%fill )  THEN
     2704                         surf_usm_v(l)%aldir(m,ind_veg_wall) =                                     &
     2705                                  building_surface_pars_f%pars(ind_s_alb_l_wall,is)
     2706                         surf_usm_v(l)%aldif(m,ind_veg_wall) =                                     &
     2707                                  building_surface_pars_f%pars(ind_s_alb_l_wall,is)
     2708                         surf_usm_v(l)%albedo_type(m,ind_veg_wall) = 0
     2709                      ENDIF
     2710
     2711                      IF ( building_surface_pars_f%pars(ind_s_alb_s_wall,is) /=                    &
     2712                           building_surface_pars_f%fill )  THEN
     2713                         surf_usm_v(l)%asdir(m,ind_veg_wall) =                                     &
     2714                                  building_surface_pars_f%pars(ind_s_alb_s_wall,is)
     2715                         surf_usm_v(l)%asdif(m,ind_veg_wall) =                                     &
     2716                                  building_surface_pars_f%pars(ind_s_alb_s_wall,is)
     2717                         surf_usm_v(l)%albedo_type(m,ind_veg_wall) = 0
     2718                      ENDIF
     2719
     2720                      IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=                     &
     2721                           building_surface_pars_f%fill )  THEN
     2722                         surf_usm_v(l)%albedo(m,ind_wat_win) =                                     &
     2723                                  building_surface_pars_f%pars(ind_s_alb_b_win,is)
     2724                         surf_usm_v(l)%albedo_type(m,ind_wat_win) = 0
     2725                      ENDIF
     2726
     2727                      IF ( building_surface_pars_f%pars(ind_s_alb_l_win,is) /=                     &
     2728                           building_surface_pars_f%fill )  THEN
     2729                         surf_usm_v(l)%aldir(m,ind_wat_win) =                                      &
     2730                                  building_surface_pars_f%pars(ind_s_alb_l_win,is)
     2731                         surf_usm_v(l)%aldif(m,ind_wat_win) =                                      &
     2732                                  building_surface_pars_f%pars(ind_s_alb_l_win,is)
     2733                         surf_usm_v(l)%albedo_type(m,ind_wat_win) = 0
     2734                      ENDIF
     2735
     2736                      IF ( building_surface_pars_f%pars(ind_s_alb_s_win,is) /=                     &
     2737                           building_surface_pars_f%fill )  THEN
     2738                         surf_usm_v(l)%asdir(m,ind_wat_win) =                                      &
     2739                                  building_surface_pars_f%pars(ind_s_alb_s_win,is)
     2740                         surf_usm_v(l)%asdif(m,ind_wat_win) =                                      &
     2741                                  building_surface_pars_f%pars(ind_s_alb_s_win,is)
     2742                         surf_usm_v(l)%albedo_type(m,ind_wat_win) = 0
     2743                      ENDIF
     2744
     2745                      IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=                   &
     2746                           building_surface_pars_f%fill )  THEN
     2747                         surf_usm_v(l)%albedo(m,ind_pav_green) =                                   &
     2748                                  building_surface_pars_f%pars(ind_s_alb_b_green,is)
     2749                         surf_usm_v(l)%albedo_type(m,ind_pav_green) = 0
     2750                      ENDIF
     2751
     2752                      IF ( building_surface_pars_f%pars(ind_s_alb_l_green,is) /=                   &
     2753                           building_surface_pars_f%fill )  THEN
     2754                         surf_usm_v(l)%aldir(m,ind_pav_green) =                                    &
     2755                                  building_surface_pars_f%pars(ind_s_alb_l_green,is)
     2756                         surf_usm_v(l)%aldif(m,ind_pav_green) =                                    &
     2757                                  building_surface_pars_f%pars(ind_s_alb_l_green,is)
     2758                         surf_usm_v(l)%albedo_type(m,ind_pav_green) = 0
     2759                      ENDIF
     2760
     2761                      IF ( building_surface_pars_f%pars(ind_s_alb_s_green,is) /=                   &
     2762                           building_surface_pars_f%fill )  THEN
     2763                         surf_usm_v(l)%asdir(m,ind_pav_green) =                                    &
     2764                                  building_surface_pars_f%pars(ind_s_alb_s_green,is)
     2765                         surf_usm_v(l)%asdif(m,ind_pav_green) =                                    &
     2766                                  building_surface_pars_f%pars(ind_s_alb_s_green,is)
     2767                         surf_usm_v(l)%albedo_type(m,ind_pav_green) = 0
     2768                      ENDIF
     2769
     2770                      EXIT ! Surface was found and processed
    20062771                   ENDIF
    20072772                ENDDO
    20082773             ENDDO
    2009 
    2010           ENDIF
    2011 !
    2012 !--       Read explicit albedo values from building surface pars. If present,
    2013 !--       they override all less specific albedo values and force a albedo_type
    2014 !--       to zero in order to take effect.
    2015           IF ( building_surface_pars_f%from_file )  THEN
    2016              DO  m = 1, surf_usm_h%ns
    2017                 i = surf_usm_h%i(m)
    2018                 j = surf_usm_h%j(m)
    2019                 k = surf_usm_h%k(m)
    2020 !
    2021 !--             Iterate over surfaces in column, check height and orientation
    2022                 DO  is = building_surface_pars_f%index_ji(1,j,i), &
    2023                          building_surface_pars_f%index_ji(2,j,i)
    2024                    IF ( building_surface_pars_f%coords(4,is) == -surf_usm_h%koff .AND. &
    2025                         building_surface_pars_f%coords(1,is) == k )  THEN
    2026 
    2027                       IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=      &
    2028                            building_surface_pars_f%fill )  THEN
    2029                          surf_usm_h%albedo(m,ind_veg_wall) =                         &
    2030                                   building_surface_pars_f%pars(ind_s_alb_b_wall,is)
    2031                          surf_usm_h%albedo_type(m,ind_veg_wall) = 0
    2032                       ENDIF
    2033 
    2034                       IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=       &
    2035                            building_surface_pars_f%fill )  THEN
    2036                          surf_usm_h%albedo(m,ind_wat_win) =                          &
    2037                                   building_surface_pars_f%pars(ind_s_alb_b_win,is)
    2038                          surf_usm_h%albedo_type(m,ind_wat_win) = 0
    2039                       ENDIF
    2040 
    2041                       IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=     &
    2042                            building_surface_pars_f%fill )  THEN
    2043                          surf_usm_h%albedo(m,ind_pav_green) =                        &
    2044                                   building_surface_pars_f%pars(ind_s_alb_b_green,is)
    2045                          surf_usm_h%albedo_type(m,ind_pav_green) = 0
    2046                       ENDIF
    2047 
    2048                       EXIT ! surface was found and processed
    2049                    ENDIF
    2050                 ENDDO
    2051              ENDDO
    2052 
    2053              DO  l = 0, 3
    2054                 DO  m = 1, surf_usm_v(l)%ns
    2055                    i = surf_usm_v(l)%i(m)
    2056                    j = surf_usm_v(l)%j(m)
    2057                    k = surf_usm_v(l)%k(m)
    2058 !
    2059 !--                Iterate over surfaces in column, check height and orientation
    2060                    DO  is = building_surface_pars_f%index_ji(1,j,i), &
    2061                             building_surface_pars_f%index_ji(2,j,i)
    2062                       IF ( building_surface_pars_f%coords(5,is) == -surf_usm_v(l)%joff .AND. &
    2063                            building_surface_pars_f%coords(6,is) == -surf_usm_v(l)%ioff .AND. &
    2064                            building_surface_pars_f%coords(1,is) == k )  THEN
    2065 
    2066                          IF ( building_surface_pars_f%pars(ind_s_alb_b_wall,is) /=      &
    2067                               building_surface_pars_f%fill )  THEN
    2068                             surf_usm_v(l)%albedo(m,ind_veg_wall) =                      &
    2069                                      building_surface_pars_f%pars(ind_s_alb_b_wall,is)
    2070                             surf_usm_v(l)%albedo_type(m,ind_veg_wall) = 0
    2071                          ENDIF
    2072 
    2073                          IF ( building_surface_pars_f%pars(ind_s_alb_b_win,is) /=       &
    2074                               building_surface_pars_f%fill )  THEN
    2075                             surf_usm_v(l)%albedo(m,ind_wat_win) =                       &
    2076                                      building_surface_pars_f%pars(ind_s_alb_b_win,is)
    2077                             surf_usm_v(l)%albedo_type(m,ind_wat_win) = 0
    2078                          ENDIF
    2079 
    2080                          IF ( building_surface_pars_f%pars(ind_s_alb_b_green,is) /=     &
    2081                               building_surface_pars_f%fill )  THEN
    2082                             surf_usm_v(l)%albedo(m,ind_pav_green) =                     &
    2083                                      building_surface_pars_f%pars(ind_s_alb_b_green,is)
    2084                             surf_usm_v(l)%albedo_type(m,ind_pav_green) = 0
    2085                          ENDIF
    2086 
    2087                          EXIT ! surface was found and processed
    2088                       ENDIF
    2089                    ENDDO
    2090                 ENDDO
    2091              ENDDO
    2092           ENDIF
    2093 !
    2094 !--    Initialization actions for RRTMG
    2095        ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
    2096 #if defined ( __rrtmg )
    2097 !
    2098 !--       Allocate albedos for short/longwave radiation, horizontal surfaces
    2099 !--       for wall/green/window (USM) or vegetation/pavement/water surfaces
    2100 !--       (LSM).
    2101           ALLOCATE ( surf_lsm_h%aldif(1:surf_lsm_h%ns,0:2)       )
    2102           ALLOCATE ( surf_lsm_h%aldir(1:surf_lsm_h%ns,0:2)       )
    2103           ALLOCATE ( surf_lsm_h%asdif(1:surf_lsm_h%ns,0:2)       )
    2104           ALLOCATE ( surf_lsm_h%asdir(1:surf_lsm_h%ns,0:2)       )
    2105           ALLOCATE ( surf_lsm_h%rrtm_aldif(1:surf_lsm_h%ns,0:2)  )
    2106           ALLOCATE ( surf_lsm_h%rrtm_aldir(1:surf_lsm_h%ns,0:2)  )
    2107           ALLOCATE ( surf_lsm_h%rrtm_asdif(1:surf_lsm_h%ns,0:2)  )
    2108           ALLOCATE ( surf_lsm_h%rrtm_asdir(1:surf_lsm_h%ns,0:2)  )
    2109 
    2110           ALLOCATE ( surf_usm_h%aldif(1:surf_usm_h%ns,0:2)       )
    2111           ALLOCATE ( surf_usm_h%aldir(1:surf_usm_h%ns,0:2)       )
    2112           ALLOCATE ( surf_usm_h%asdif(1:surf_usm_h%ns,0:2)       )
    2113           ALLOCATE ( surf_usm_h%asdir(1:surf_usm_h%ns,0:2)       )
    2114           ALLOCATE ( surf_usm_h%rrtm_aldif(1:surf_usm_h%ns,0:2)  )
    2115           ALLOCATE ( surf_usm_h%rrtm_aldir(1:surf_usm_h%ns,0:2)  )
    2116           ALLOCATE ( surf_usm_h%rrtm_asdif(1:surf_usm_h%ns,0:2)  )
    2117           ALLOCATE ( surf_usm_h%rrtm_asdir(1:surf_usm_h%ns,0:2)  )
    2118 
    2119 !
    2120 !--       Allocate broadband albedo (temporary for the current radiation
    2121 !--       implementations)
    2122           IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )                            &
    2123              ALLOCATE( surf_lsm_h%albedo(1:surf_lsm_h%ns,0:2)     )
    2124           IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )                            &
    2125              ALLOCATE( surf_usm_h%albedo(1:surf_usm_h%ns,0:2)     )
    2126 
    2127 !
    2128 !--       Allocate albedos for short/longwave radiation, vertical surfaces
     2774          ENDDO
     2775       ENDIF
     2776
     2777!
     2778!--    Calculate initial values of current (cosine of) zenith angle and whether the sun is up
     2779       CALL get_date_time( time_since_reference_point, day_of_year = day_of_year,                  &
     2780                           second_of_day = second_of_day )
     2781       CALL calc_zenith( day_of_year, second_of_day )
     2782!
     2783!--    Calculate initial surface albedo for different surfaces
     2784       IF ( .NOT. constant_albedo )  THEN
     2785#if defined( __netcdf )
     2786!
     2787!--       Horizontally aligned natural and urban surfaces
     2788          CALL calc_albedo( surf_lsm_h )
     2789          CALL calc_albedo( surf_usm_h )
     2790!
     2791!--       Vertically aligned natural and urban surfaces
    21292792          DO  l = 0, 3
    2130 
    2131              ALLOCATE ( surf_lsm_v(l)%aldif(1:surf_lsm_v(l)%ns,0:2)      )
    2132              ALLOCATE ( surf_lsm_v(l)%aldir(1:surf_lsm_v(l)%ns,0:2)      )
    2133              ALLOCATE ( surf_lsm_v(l)%asdif(1:surf_lsm_v(l)%ns,0:2)      )
    2134              ALLOCATE ( surf_lsm_v(l)%asdir(1:surf_lsm_v(l)%ns,0:2)      )
    2135 
    2136              ALLOCATE ( surf_lsm_v(l)%rrtm_aldif(1:surf_lsm_v(l)%ns,0:2) )
    2137              ALLOCATE ( surf_lsm_v(l)%rrtm_aldir(1:surf_lsm_v(l)%ns,0:2) )
    2138              ALLOCATE ( surf_lsm_v(l)%rrtm_asdif(1:surf_lsm_v(l)%ns,0:2) )
    2139              ALLOCATE ( surf_lsm_v(l)%rrtm_asdir(1:surf_lsm_v(l)%ns,0:2) )
    2140 
    2141              ALLOCATE ( surf_usm_v(l)%aldif(1:surf_usm_v(l)%ns,0:2)      )
    2142              ALLOCATE ( surf_usm_v(l)%aldir(1:surf_usm_v(l)%ns,0:2)      )
    2143              ALLOCATE ( surf_usm_v(l)%asdif(1:surf_usm_v(l)%ns,0:2)      )
    2144              ALLOCATE ( surf_usm_v(l)%asdir(1:surf_usm_v(l)%ns,0:2)      )
    2145 
    2146              ALLOCATE ( surf_usm_v(l)%rrtm_aldif(1:surf_usm_v(l)%ns,0:2) )
    2147              ALLOCATE ( surf_usm_v(l)%rrtm_aldir(1:surf_usm_v(l)%ns,0:2) )
    2148              ALLOCATE ( surf_usm_v(l)%rrtm_asdif(1:surf_usm_v(l)%ns,0:2) )
    2149              ALLOCATE ( surf_usm_v(l)%rrtm_asdir(1:surf_usm_v(l)%ns,0:2) )
    2150 !
    2151 !--          Allocate broadband albedo (temporary for the current radiation
    2152 !--          implementations)
    2153              IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )                    &
    2154                 ALLOCATE( surf_lsm_v(l)%albedo(1:surf_lsm_v(l)%ns,0:2) )
    2155              IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )                    &
    2156                 ALLOCATE( surf_usm_v(l)%albedo(1:surf_usm_v(l)%ns,0:2) )
    2157 
     2793             CALL calc_albedo( surf_lsm_v(l) )
     2794             CALL calc_albedo( surf_usm_v(l) )
    21582795          ENDDO
    2159 !
    2160 !--       Level 1 initialization of spectral albedos via namelist
    2161 !--       paramters. Please note, this case all surface tiles are initialized
    2162 !--       the same.
     2796#endif
     2797       ELSE
     2798!
     2799!--       Initialize sun-inclination independent spectral albedos
     2800!--       Horizontal surfaces
    21632801          IF ( surf_lsm_h%ns > 0 )  THEN
    2164              surf_lsm_h%aldif  = albedo_lw_dif
    2165              surf_lsm_h%aldir  = albedo_lw_dir
    2166              surf_lsm_h%asdif  = albedo_sw_dif
    2167              surf_lsm_h%asdir  = albedo_sw_dir