Changeset 4634 for palm


Ignore:
Timestamp:
Aug 5, 2020 3:42:28 PM (3 months ago)
Author:
pavelkrc
Message:

Revert re-formatting (r4624) which cannot be merged with current development (see next revisions)

File:
1 edited

Legend:

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

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