Changeset 4452

Timestamp:

Mar 10, 2020 8:15:32 PM (5 years ago)

Author:

suehring

Message:

Bugfix in calc_albedo

File:

• palm/trunk/SOURCE/radiation_model_mod.f90 (modified) (334 diffs)

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $Id$
-! - Change order of dimension in surface arrays %frac, %emissivity and %albedo 
+! Bugfix in calc_albedo
+!
+! 4442 2020-03-04 19:21:13Z suehring
+! - Change order of dimension in surface arrays %frac, %emissivity and %albedo
 !   to allow for better vectorization in the radiation interactions.
 ! - Minor formatting issues
-! 
+!
 ! 4441 2020-03-04 19:20:35Z suehring
 ! bugfixes: cpp-directives for serial mode moved, small changes to get serial mode compiled
-! 
+!
 ! 4400 2020-02-10 20:32:41Z suehring
 ! Initialize radiation arrays with zero
-! 
+!
 ! 4392 2020-01-31 16:14:57Z pavelkrc
 ! - Add debug tracing of large radiative fluxes (option trace_fluxes_above)
 ! - Print exact counts of SVF and CSF if debut_output is enabled
 ! - Update descriptions of RTM 3.0 and related comments
-! 
+!
 ! 4360 2020-01-07 11:25:50Z suehring
 ! Renamed pc_heating_rate, pc_transpiration_rate, pc_transpiration_rate to
 ! pcm_heating_rate, pcm_latent_rate, pcm_transpiration_rate
-! 
+!
 ! 4340 2019-12-16 08:17:03Z Giersch
-! Albedo indices for building_surface_pars are now declared as parameters to 
+! Albedo indices for building_surface_pars are now declared as parameters to
 ! prevent an error if the gfortran compiler with -Werror=unused-value is used
-!  
+!
 ! 4291 2019-11-11 12:36:54Z moh.hefny
-! Enabled RTM in case of biometeorology even if there is no vertical 
+! Enabled RTM in case of biometeorology even if there is no vertical
 ! surfaces or 3D vegetation in the domain
-! 
+!
 ! 4286 2019-10-30 16:01:14Z resler
 ! - Fix wrong treating of time_rad during interpolation in radiation_model_mod
@@ -62 +65 @@
 ! 4271 2019-10-23 10:46:41Z maronga
 ! Bugfix: missing parentheses in calculation of snow albedo
-! 
+!
 ! 4245 2019-09-30 08:40:37Z pavelkrc
 ! Initialize explicit per-surface albedos from building_surface_pars
-! 
+!
 ! 4238 2019-09-25 16:06:01Z suehring
 ! Modify check in order to avoid equality comparisons of floating points
-! 
+!
 ! 4227 2019-09-10 18:04:34Z gronemeier
 ! implement new palm_date_time_mod
-! 
+!
 ! 4226 2019-09-10 17:03:24Z suehring
 ! - Netcdf input routine for dimension length renamed
-! - Define time variable for external radiation input relative to time_utc_init 
-! 
+! - Define time variable for external radiation input relative to time_utc_init
+!
 ! 4210 2019-09-02 13:07:09Z suehring
 ! - Revise steering of splitting diffuse and direct radiation
@@ -81 +84 @@
 ! - Optimize mapping of radiation components onto 2D arrays, avoid unnecessary
 !   operations
-! 
+!
 ! 4208 2019-09-02 09:01:07Z suehring
 ! Bugfix in accessing albedo_pars in the clear-sky branch
-! (merge from branch resler) 
-! 
+! (merge from branch resler)
+!
 ! 4198 2019-08-29 15:17:48Z gronemeier
 ! Prohibit execution of radiation model if rotation_angle is not zero
-! 
+!
 ! 4197 2019-08-29 14:33:32Z suehring
 ! Revise steering of surface albedo initialization when albedo_pars is provided
-! 
+!
 ! 4190 2019-08-27 15:42:37Z suehring
-! Implement external radiation forcing also for level-of-detail = 2 
+! Implement external radiation forcing also for level-of-detail = 2
 ! (horizontally 2D radiation)
-! 
+!
 ! 4188 2019-08-26 14:15:47Z suehring
 ! Minor adjustment in error message
-! 
+!
 ! 4187 2019-08-26 12:43:15Z suehring
-! - Take external radiation from root domain dynamic input if not provided for 
+! - Take external radiation from root domain dynamic input if not provided for
 !   each nested domain
 ! - Combine MPI_ALLREDUCE calls to reduce mpi overhead
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 4179 2019-08-21 11:16:12Z suehring
 ! Remove debug prints
-! 
+!
 ! 4178 2019-08-21 11:13:06Z suehring
 ! External radiation forcing implemented.
-! 
+!
 ! 4168 2019-08-16 13:50:17Z suehring
 ! Replace function get_topography_top_index by topo_top_ind
-! 
+!
 ! 4157 2019-08-14 09:19:12Z suehring
 ! Give informative message on raytracing distance only by core zero
-! 
+!
 ! 4148 2019-08-08 11:26:00Z suehring
 ! Comments added
-! 
+!
 ! 4134 2019-08-02 18:39:57Z suehring
 ! Bugfix in formatted write statement
-! 
+!
 ! 4127 2019-07-30 14:47:10Z suehring
 ! Remove unused pch_index (merge from branch resler)
-! 
+!
 ! 4089 2019-07-11 14:30:27Z suehring
 ! - Correct level 2 initialization of spectral albedos in rrtmg branch, long- and
-!   shortwave albedos were mixed-up. 
+!   shortwave albedos were mixed-up.
 ! - Change order of albedo_pars so that it is now consistent with the defined
-!   order of albedo_pars in PIDS 
-! 
+!   order of albedo_pars in PIDS
+!
 ! 4069 2019-07-01 14:05:51Z Giersch
-! Masked output running index mid has been introduced as a local variable to 
+! Masked output running index mid has been introduced as a local variable to
 ! avoid runtime error (Loop variable has been modified) in time_integration
-! 
+!
 ! 4067 2019-07-01 13:29:25Z suehring
 ! Bugfix, pass dummy string to MPI_INFO_SET (J. Resler)
-! 
+!
 ! 4039 2019-06-18 10:32:41Z suehring
 ! Bugfix for masked data output
-! 
+!
 ! 4008 2019-05-30 09:50:11Z moh.hefny
 ! Bugfix in check variable when a variable's string is less than 3
@@ -150 +153 @@
 !
 ! 3992 2019-05-22 16:49:38Z suehring
-! Bugfix in rrtmg radiation branch in a nested run when the lowest prognistic 
-! grid points in a child domain are all inside topography 
-! 
+! Bugfix in rrtmg radiation branch in a nested run when the lowest prognistic
+! grid points in a child domain are all inside topography
+!
 ! 3987 2019-05-22 09:52:13Z kanani
 ! Introduce alternative switch for debug output during timestepping
-! 
+!
 ! 3943 2019-05-02 09:50:41Z maronga
 ! Missing blank characteer added.
-! 
+!
 ! 3900 2019-04-16 15:17:43Z suehring
 ! Fixed initialization problem
-! 
+!
 ! 3885 2019-04-11 11:29:34Z kanani
-! Changes related to global restructuring of location messages and introduction 
+! Changes related to global restructuring of location messages and introduction
 ! of additional debug messages
-! 
+!
 ! 3881 2019-04-10 09:31:22Z suehring
 ! Output of albedo and emissivity moved from USM, bugfixes in initialization
 ! of albedo
-! 
+!
 ! 3861 2019-04-04 06:27:41Z maronga
 ! Bugfix: factor of 4.0 required instead of 3.0 in calculation of rad_lw_out_change_0
-! 
+!
 ! 3859 2019-04-03 20:30:31Z maronga
 ! Added some descriptions
-! 
+!
 ! 3847 2019-04-01 14:51:44Z suehring
 ! Implement check for dt_radiation (must be > 0)
-! 
+!
 ! 3846 2019-04-01 13:55:30Z suehring
 ! unused variable removed
-! 
+!
 ! 3814 2019-03-26 08:40:31Z pavelkrc
 ! Change zenith(0:0) and others to scalar.
 ! Code review.
 ! Rename exported nzu, nzp and related variables due to name conflict
-! 
+!
 ! 3771 2019-02-28 12:19:33Z raasch
 ! rrtmg preprocessor for directives moved/added, save attribute added to temporary
 ! pointers to avoid compiler warnings about outlived pointer targets,
 ! statement added to avoid compiler warning about unused variable
-! 
+!
 ! 3769 2019-02-28 10:16:49Z moh.hefny
 ! removed unused variables and subroutine radiation_radflux_gridbox
@@ -197 +200 @@
 ! 3767 2019-02-27 08:18:02Z raasch
 ! unused variable for file index removed from rrd-subroutines parameter list
-! 
+!
 ! 3760 2019-02-21 18:47:35Z moh.hefny
 ! Bugfix: initialized simulated_time before calculating solar position
@@ -205 +208 @@
 ! (resler, pavelkrc)
 ! Bugfixes: add further required MRT factors to read/write_svf,
-! fix for aggregating view factors to eliminate local noise in reflected 
-! irradiance at mutually close surfaces (corners, presence of trees) in the 
+! fix for aggregating view factors to eliminate local noise in reflected
+! irradiance at mutually close surfaces (corners, presence of trees) in the
 ! angular discretization scheme.
 !
@@ -213 +216 @@
 !
 ! 3705 2019-01-29 19:56:39Z suehring
-! Make variables that are sampled in virtual measurement module public 
-! 
+! Make variables that are sampled in virtual measurement module public
+!
 ! 3704 2019-01-29 19:51:41Z suehring
 ! Some interface calls moved to module_interface + cleanup
-! 
+!
 ! 3667 2019-01-10 14:26:24Z schwenkel
 ! Modified check for rrtmg input files
-! 
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! nopointer option removed
-! 
+!
 ! 1496 2014-12-02 17:25:50Z maronga
 ! Initial revision
-! 
+!
 !
 ! Description:
@@ -247 +250 @@
 !------------------------------------------------------------------------------!
  MODULE radiation_model_mod
- 
+
     USE arrays_3d,                                                             &
         ONLY:  dzw, hyp, nc, pt, p, q, ql, u, v, w, zu, zw, exner, d_exner
@@ -276 +279 @@
 
     USE grid_variables,                                                        &
-         ONLY:  ddx, ddy, dx, dy 
+         ONLY:  ddx, ddy, dx, dy
 
     USE indices,                                                               &
@@ -366 +369 @@
 !-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
     CHARACTER(37), DIMENSION(0:33), PARAMETER :: albedo_type_name = (/      &
-                                   'user defined                         ', & !  0 
+                                   'user defined                         ', & !  0
                                    'ocean                                ', & !  1
                                    'mixed farming, tall grassland        ', & !  2
-                                   'tall/medium grassland                ', & !  3 
+                                   'tall/medium grassland                ', & !  3
                                    'evergreen shrubland                  ', & !  4
                                    'short grassland/meadow/shrubland     ', & !  5
@@ -377 +380 @@
                                    'tropical evergreen broadleaved forest', & !  9
                                    'medium/tall grassland/woodland       ', & ! 10
-                                   'desert, sandy                        ', & ! 11 
-                                   'desert, rocky                        ', & ! 12 
+                                   'desert, sandy                        ', & ! 11
+                                   'desert, rocky                        ', & ! 12
                                    'tundra                               ', & ! 13
-                                   'land ice                             ', & ! 14 
-                                   'sea ice                              ', & ! 15 
+                                   'land ice                             ', & ! 14
+                                   'sea ice                              ', & ! 15
                                    'snow                                 ', & ! 16
                                    'bare soil                            ', & ! 17
@@ -415 +418 @@
 
     INTEGER(iwp) :: albedo_type  = 9999999_iwp, &     !< Albedo surface type
-                    dots_rad     = 0_iwp              !< starting index for timeseries output 
+                    dots_rad     = 0_iwp              !< starting index for timeseries output
 
     LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
@@ -467 +470 @@
     REAL(wp), PARAMETER :: emissivity_atm_clsky = 0.8_wp       !< emissivity of the clear-sky atmosphere
 !
-!-- Land surface albedos for solar zenith angle of 60degree after Briegleb (1992)     
+!-- Land surface albedos for solar zenith angle of 60degree after Briegleb (1992)
 !-- (broadband, longwave, shortwave ):   bb,      lw,      sw,
-    REAL(wp), DIMENSION(0:2,1:33), PARAMETER :: albedo_pars = RESHAPE( (/& 
+    REAL(wp), DIMENSION(0:2,1:33), PARAMETER :: albedo_pars = RESHAPE( (/&
                                    0.06_wp, 0.06_wp, 0.06_wp,            & !  1
                                    0.19_wp, 0.28_wp, 0.09_wp,            & !  2
@@ -541 +544 @@
 
 !
-!-- The following variables should be only changed with care, as this will 
+!-- The following variables should be only changed with care, as this will
 !-- require further setting of some variables, which is currently not
 !-- implemented (aerosols, ice phase).
@@ -570 +573 @@
                                              rrtm_cliqwp,    & !< in-cloud liquid water path (g/m2)
                                              rrtm_co2vmr,    & !< CO2 volume mixing ratio (g/mol)
-                                             rrtm_emis,      & !< surface emissivity (0-1)  
+                                             rrtm_emis,      & !< surface emissivity (0-1)
                                              rrtm_h2ovmr,    & !< H2O volume mixing ratio
                                              rrtm_n2ovmr,    & !< N2O volume mixing ratio
@@ -582 +585 @@
                                              rrtm_tlev,      & !< actual temperature (K, zw-grid)
                                              rrtm_lwdflx,    & !< RRTM output of incoming longwave radiation flux (W/m2)
-                                             rrtm_lwdflxc,   & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2) 
+                                             rrtm_lwdflxc,   & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
                                              rrtm_lwuflx,    & !< RRTM output of outgoing longwave radiation flux (W/m2)
                                              rrtm_lwuflxc,   & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
@@ -590 +593 @@
                                              rrtm_lwhrc,     & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
                                              rrtm_swdflx,    & !< RRTM output of incoming shortwave radiation flux (W/m2)
-                                             rrtm_swdflxc,   & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2) 
+                                             rrtm_swdflxc,   & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
                                              rrtm_swuflx,    & !< RRTM output of outgoing shortwave radiation flux (W/m2)
                                              rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
                                              rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
-                                             rrtm_swhrc,     & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d) 
-                                             rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m2) 
-                                             rrtm_difdflux     !< RRTM output of incoming diffuse shortwave (W/m2) 
+                                             rrtm_swhrc,     & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
+                                             rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m2)
+                                             rrtm_difdflux     !< RRTM output of incoming diffuse shortwave (W/m2)
 
     REAL(wp), DIMENSION(1) ::                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
@@ -734 +737 @@
 !-- type for calculation of svf
     TYPE t_svf
-        INTEGER(iwp)                               :: isurflt           !< 
-        INTEGER(iwp)                               :: isurfs            !< 
-        REAL(wp)                                   :: rsvf              !< 
-        REAL(wp)                                   :: rtransp           !< 
+        INTEGER(iwp)                               :: isurflt           !<
+        INTEGER(iwp)                               :: isurfs            !<
+        REAL(wp)                                   :: rsvf              !<
+        REAL(wp)                                   :: rtransp           !<
     END TYPE
 
 !-- type for calculation of csf
     TYPE t_csf
-        INTEGER(iwp)                               :: ip                !< 
-        INTEGER(iwp)                               :: itx               !< 
-        INTEGER(iwp)                               :: ity               !< 
-        INTEGER(iwp)                               :: itz               !< 
+        INTEGER(iwp)                               :: ip                !<
+        INTEGER(iwp)                               :: itx               !<
+        INTEGER(iwp)                               :: ity               !<
+        INTEGER(iwp)                               :: itz               !<
         INTEGER(iwp)                               :: isurfs            !< Idx of source face / -1 for sky
         REAL(wp)                                   :: rcvf              !< Canopy view factor for faces /
@@ -788 +791 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif      !< array of diffuse sw radiation from sky and model boundary falling to local surface
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif      !< array of diffuse lw radiation from sky and model boundary falling to local surface
-    
+
                                                                         !< Outward radiation is only valid for nonvirtual surfaces
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsl        !< array of reflected sw radiation for local surface in i-th reflection
@@ -835 +838 @@
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  boxes            !< coordinates of gridboxes being crossed by ray
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  crlens           !< array of crossing lengths of ray for particular grid boxes
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  lad_ip           !< array of numbers of process where lad is stored 
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  lad_ip           !< array of numbers of process where lad is stored
 #if defined( __parallel )
     INTEGER(kind=MPI_ADDRESS_KIND), &
@@ -876 +879 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE            :: emiss_surf         !< emissivity of the wall surface
 !
-!-- External radiation. Depending on the given level of detail either a 1D or 
+!-- External radiation. Depending on the given level of detail either a 1D or
 !-- a 3D array will be allocated.
-    TYPE( real_1d_3d ) ::  rad_lw_in_f     !< external incoming longwave radiation, from observation or model 
-    TYPE( real_1d_3d ) ::  rad_sw_in_f     !< external incoming shortwave radiation, from observation or model 
+    TYPE( real_1d_3d ) ::  rad_lw_in_f     !< external incoming longwave radiation, from observation or model
+    TYPE( real_1d_3d ) ::  rad_sw_in_f     !< external incoming shortwave radiation, from observation or model
     TYPE( real_1d_3d ) ::  rad_sw_in_dif_f !< external incoming shortwave radiation, diffuse part, from observation or model
-    TYPE( real_1d_3d ) ::  time_rad_f      !< time dimension for external radiation, from observation or model 
+    TYPE( real_1d_3d ) ::  time_rad_f      !< time dimension for external radiation, from observation or model
 
     INTERFACE radiation_check_data_output
@@ -894 +897 @@
        MODULE PROCEDURE radiation_check_data_output_pr
     END INTERFACE radiation_check_data_output_pr
-  
+
     INTERFACE radiation_check_parameters
        MODULE PROCEDURE radiation_check_parameters
     END INTERFACE radiation_check_parameters
-  
+
     INTERFACE radiation_clearsky
        MODULE PROCEDURE radiation_clearsky
     END INTERFACE radiation_clearsky
- 
+
     INTERFACE radiation_constant
        MODULE PROCEDURE radiation_constant
     END INTERFACE radiation_constant
- 
+
     INTERFACE radiation_control
        MODULE PROCEDURE radiation_control
@@ -933 +936 @@
     INTERFACE radiation_header
        MODULE PROCEDURE radiation_header
-    END INTERFACE radiation_header 
-  
+    END INTERFACE radiation_header
+
     INTERFACE radiation_init
        MODULE PROCEDURE radiation_init
@@ -942 +945 @@
        MODULE PROCEDURE radiation_parin
     END INTERFACE radiation_parin
-    
+
     INTERFACE radiation_rrtmg
        MODULE PROCEDURE radiation_rrtmg
@@ -969 +972 @@
        MODULE PROCEDURE radiation_interaction_init
     END INTERFACE radiation_interaction_init
- 
+
     INTERFACE radiation_presimulate_solar_pos
        MODULE PROCEDURE radiation_presimulate_solar_pos
@@ -1005 +1008 @@
            radiation_read_svf, radiation_presimulate_solar_pos
 
-    
+
 !
 !-- Public variables and constants / NEEDS SORTING
@@ -1042 +1045 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE radiation_control
- 
- 
+
+
        IMPLICIT NONE
 
@@ -1054 +1057 @@
           CASE ( 'constant' )
              CALL radiation_constant
-          
-          CASE ( 'clear-sky' ) 
+
+          CASE ( 'clear-sky' )
              CALL radiation_clearsky
-       
+
           CASE ( 'rrtmg' )
              CALL radiation_rrtmg
-             
+
           CASE ( 'external' )
 !
@@ -1084 +1087 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE radiation_check_data_output( variable, unit, i, ilen, k )
- 
- 
+
+
        USE control_parameters,                                                 &
            ONLY: data_output, message_string
@@ -1262 +1265 @@
 ! ------------
 !> Check data output of profiles for radiation model
-!------------------------------------------------------------------------------!  
+!------------------------------------------------------------------------------!
     SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit,      &
                dopr_unit )
- 
+
        USE arrays_3d,                                                          &
            ONLY: zu
@@ -1279 +1282 @@
 
        IMPLICIT NONE
-   
-       CHARACTER (LEN=*) ::  unit      !< 
-       CHARACTER (LEN=*) ::  variable  !< 
+
+       CHARACTER (LEN=*) ::  unit      !<
+       CHARACTER (LEN=*) ::  variable  !<
        CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
- 
-       INTEGER(iwp) ::  var_count     !< 
+
+       INTEGER(iwp) ::  var_count     !<
 
        SELECT CASE ( TRIM( variable ) )
-       
+
          CASE ( 'rad_net' )
              IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
@@ -1315 +1318 @@
                 dopr_unit  = 'W/m2'
                 hom(:,2,100,:)  = SPREAD( zw, 2, statistic_regions+1 )
-                unit = dopr_unit 
+                unit = dopr_unit
              ENDIF
 
@@ -1330 +1333 @@
                 dopr_unit  = 'W/m2'
                 hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions+1 )
-                unit = dopr_unit   
+                unit = dopr_unit
              ENDIF
 
@@ -1431 +1434 @@
 
     END SUBROUTINE radiation_check_data_output_pr
- 
- 
+
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -1444 +1447 @@
 
        USE netcdf_data_input_mod,                                              &
-           ONLY:  input_pids_static                 
-    
+           ONLY:  input_pids_static
+
        IMPLICIT NONE
-       
-!
-!--    In case no urban-surface or land-surface model is applied, usage of 
-!--    a radiation model make no sense.         
+
+!
+!--    In case no urban-surface or land-surface model is applied, usage of
+!--    a radiation model make no sense.
        IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
           message_string = 'Usage of radiation module is only allowed if ' //  &
@@ -1466 +1469 @@
        ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
 #if ! defined ( __rrtmg )
-          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
+          message_string = 'radiation_scheme = "rrtmg" requires ' //           &
                            'compilation of PALM with pre-processor ' //        &
                            'directive -D__rrtmg'
@@ -1472 +1475 @@
 #endif
 #if defined ( __rrtmg ) && ! defined( __netcdf )
-          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
+          message_string = 'radiation_scheme = "rrtmg" requires ' //           &
                            'the use of NetCDF (preprocessor directive ' //     &
                            '-D__netcdf'
@@ -1480 +1483 @@
        ENDIF
 !
-!--    Checks performed only if data is given via namelist only. 
+!--    Checks performed only if data is given via namelist only.
        IF ( .NOT. input_pids_static )  THEN
           IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.          &
                radiation_scheme == 'clear-sky')  THEN
-             message_string = 'radiation_scheme = "clear-sky" in combination'//& 
+             message_string = 'radiation_scheme = "clear-sky" in combination'//&
                               'with albedo_type = 0 requires setting of'//     &
                               'albedo /= 9999999.9'
@@ -1492 +1495 @@
           IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.     &
              ( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
-          .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& 
+          .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp&
              ) ) THEN
-             message_string = 'radiation_scheme = "rrtmg" in combination' //   & 
+             message_string = 'radiation_scheme = "rrtmg" in combination' //   &
                               'with albedo_type = 0 requires setting of ' //   &
                               'albedo_lw_dif /= 9999999.9' //                  &
@@ -1506 +1509 @@
 !--    Parallel rad_angular_discretization without raytrace_mpi_rma is not implemented
 !--    Serial mode does not allow mpi_rma
-#if defined( __parallel )      
+#if defined( __parallel )
        IF ( rad_angular_discretization  .AND.  .NOT. raytrace_mpi_rma )  THEN
           message_string = 'rad_angular_discretization can only be used ' //  &
@@ -1522 +1525 @@
        IF ( cloud_droplets  .AND.   radiation_scheme == 'rrtmg'  .AND.         &
             average_radiation ) THEN
-          message_string = 'average_radiation = .T. with radiation_scheme'//   & 
+          message_string = 'average_radiation = .T. with radiation_scheme'//   &
                            '= "rrtmg" in combination cloud_droplets = .T.'//   &
                            'is not implementd'
@@ -1539 +1542 @@
 !--    Check for dt_radiation
        IF ( dt_radiation <= 0.0 )  THEN
-          message_string = 'dt_radiation must be > 0.0'  
-          CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 ) 
+          message_string = 'dt_radiation must be > 0.0'
+          CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 )
        ENDIF
 !
@@ -1549 +1552 @@
                            'model.&Using rotation_angle /= 0.0 is not allowed in combination ' //  &
                            'with the radiation model at the moment!'
-          CALL message( 'check_parameters', 'PA0675', 1, 2, 0, 6, 0 ) 
+          CALL message( 'check_parameters', 'PA0675', 1, 2, 0, 6, 0 )
        ENDIF
- 
-    END SUBROUTINE radiation_check_parameters 
- 
- 
+
+    END SUBROUTINE radiation_check_parameters
+
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -1561 +1564 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE radiation_init
-    
+
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i         !< running index x-direction 
+       INTEGER(iwp) ::  i         !< running index x-direction
        INTEGER(iwp) ::  is        !< running index for input surface elements
-       INTEGER(iwp) ::  ioff      !< offset in x between surface element reference grid point in atmosphere and actual surface 
-       INTEGER(iwp) ::  j         !< running index y-direction 
-       INTEGER(iwp) ::  joff      !< offset in y between surface element reference grid point in atmosphere and actual surface 
+       INTEGER(iwp) ::  ioff      !< offset in x between surface element reference grid point in atmosphere and actual surface
+       INTEGER(iwp) ::  j         !< running index y-direction
+       INTEGER(iwp) ::  joff      !< offset in y between surface element reference grid point in atmosphere and actual surface
        INTEGER(iwp) ::  k         !< running index z-direction
-       INTEGER(iwp) ::  l         !< running index for orientation of vertical surfaces 
+       INTEGER(iwp) ::  l         !< running index for orientation of vertical surfaces
        INTEGER(iwp) ::  m         !< running index for surface elements
        INTEGER(iwp) ::  ntime = 0 !< number of available external radiation timesteps
@@ -1584 +1587 @@
 !      or if biometeorology output is required for flat surfaces.
 !--    The namelist parameter radiation_interactions_on can override this behavior.
-!--    (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in 
+!--    (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in
 !--    init_surface_arrays.)
        IF ( radiation_interactions_on )  THEN
@@ -1591 +1594 @@
              average_radiation         = .TRUE.
           ELSE
-             radiation_interactions_on = .FALSE.   !< reset namelist parameter: no interactions 
+             radiation_interactions_on = .FALSE.   !< reset namelist parameter: no interactions
                                                    !< calculations necessary in case of flat surface
           ENDIF
        ELSEIF ( vertical_surfaces_exist  .OR.  plant_canopy  .OR.  biometeorology )  THEN
           message_string = 'radiation_interactions_on is set to .FALSE. although '     // &
-                           'vertical surfaces and/or trees or biometeorology exist '   // & 
+                           'vertical surfaces and/or trees or biometeorology exist '   // &
                            'is ON. The model will run without RTM (no shadows, no '    // &
                            'radiation reflections)'
@@ -1607 +1610 @@
 
 !
-!--    If required, initialize radiation interactions between surfaces 
+!--    If required, initialize radiation interactions between surfaces
 !--    via sky-view factors. This must be done before radiation is initialized.
        IF ( radiation_interactions )  CALL radiation_interaction_init
@@ -1615 +1618 @@
                   surf_lsm_h%ns > 0  )   THEN
           ALLOCATE( surf_lsm_h%rad_net(1:surf_lsm_h%ns) )
-          surf_lsm_h%rad_net = 0.0_wp 
+          surf_lsm_h%rad_net = 0.0_wp
        ENDIF
        IF ( .NOT. ALLOCATED ( surf_usm_h%rad_net )  .AND.                      &
                   surf_usm_h%ns > 0  )  THEN
           ALLOCATE( surf_usm_h%rad_net(1:surf_usm_h%ns) )
-          surf_usm_h%rad_net = 0.0_wp 
+          surf_usm_h%rad_net = 0.0_wp
        ENDIF
        DO  l = 0, 3
@@ -1626 +1629 @@
                      surf_lsm_v(l)%ns > 0  )  THEN
              ALLOCATE( surf_lsm_v(l)%rad_net(1:surf_lsm_v(l)%ns) )
-             surf_lsm_v(l)%rad_net = 0.0_wp 
+             surf_lsm_v(l)%rad_net = 0.0_wp
           ENDIF
           IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_net )  .AND.                &
                      surf_usm_v(l)%ns > 0  )  THEN
              ALLOCATE( surf_usm_v(l)%rad_net(1:surf_usm_v(l)%ns) )
-             surf_usm_v(l)%rad_net = 0.0_wp 
+             surf_usm_v(l)%rad_net = 0.0_wp
           ENDIF
        ENDDO
@@ -1641 +1644 @@
                   surf_lsm_h%ns > 0  )   THEN
           ALLOCATE( surf_lsm_h%rad_lw_out_change_0(1:surf_lsm_h%ns) )
-          surf_lsm_h%rad_lw_out_change_0 = 0.0_wp 
+          surf_lsm_h%rad_lw_out_change_0 = 0.0_wp
        ENDIF
        IF ( .NOT. ALLOCATED ( surf_usm_h%rad_lw_out_change_0 )  .AND.          &
                   surf_usm_h%ns > 0  )  THEN
           ALLOCATE( surf_usm_h%rad_lw_out_change_0(1:surf_usm_h%ns) )
-          surf_usm_h%rad_lw_out_change_0 = 0.0_wp 
+          surf_usm_h%rad_lw_out_change_0 = 0.0_wp
        ENDIF
        DO  l = 0, 3
@@ -1652 +1655 @@
                      surf_lsm_v(l)%ns > 0  )  THEN
              ALLOCATE( surf_lsm_v(l)%rad_lw_out_change_0(1:surf_lsm_v(l)%ns) )
-             surf_lsm_v(l)%rad_lw_out_change_0 = 0.0_wp 
+             surf_lsm_v(l)%rad_lw_out_change_0 = 0.0_wp
           ENDIF
           IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_lw_out_change_0 )  .AND.    &
                      surf_usm_v(l)%ns > 0  )  THEN
              ALLOCATE( surf_usm_v(l)%rad_lw_out_change_0(1:surf_usm_v(l)%ns) )
-             surf_usm_v(l)%rad_lw_out_change_0 = 0.0_wp 
+             surf_usm_v(l)%rad_lw_out_change_0 = 0.0_wp
           ENDIF
        ENDDO
@@ -1676 +1679 @@
           ALLOCATE( surf_lsm_h%rad_lw_ref(1:surf_lsm_h%ns) )
           ALLOCATE( surf_lsm_h%rad_lw_res(1:surf_lsm_h%ns) )
-          surf_lsm_h%rad_sw_in  = 0.0_wp 
-          surf_lsm_h%rad_sw_out = 0.0_wp 
-          surf_lsm_h%rad_sw_dir = 0.0_wp 
-          surf_lsm_h%rad_sw_dif = 0.0_wp 
-          surf_lsm_h%rad_sw_ref = 0.0_wp 
-          surf_lsm_h%rad_sw_res = 0.0_wp 
-          surf_lsm_h%rad_lw_in  = 0.0_wp 
-          surf_lsm_h%rad_lw_out = 0.0_wp 
-          surf_lsm_h%rad_lw_dif = 0.0_wp 
-          surf_lsm_h%rad_lw_ref = 0.0_wp 
-          surf_lsm_h%rad_lw_res = 0.0_wp 
+          surf_lsm_h%rad_sw_in  = 0.0_wp
+          surf_lsm_h%rad_sw_out = 0.0_wp
+          surf_lsm_h%rad_sw_dir = 0.0_wp
+          surf_lsm_h%rad_sw_dif = 0.0_wp
+          surf_lsm_h%rad_sw_ref = 0.0_wp
+          surf_lsm_h%rad_sw_res = 0.0_wp
+          surf_lsm_h%rad_lw_in  = 0.0_wp
+          surf_lsm_h%rad_lw_out = 0.0_wp
+          surf_lsm_h%rad_lw_dif = 0.0_wp
+          surf_lsm_h%rad_lw_ref = 0.0_wp
+          surf_lsm_h%rad_lw_res = 0.0_wp
        ENDIF
        IF ( .NOT. ALLOCATED ( surf_usm_h%rad_sw_in )  .AND.                    &
@@ -1701 +1704 @@
           ALLOCATE( surf_usm_h%rad_lw_ref(1:surf_usm_h%ns) )
           ALLOCATE( surf_usm_h%rad_lw_res(1:surf_usm_h%ns) )
-          surf_usm_h%rad_sw_in  = 0.0_wp 
-          surf_usm_h%rad_sw_out = 0.0_wp 
-          surf_usm_h%rad_sw_dir = 0.0_wp 
-          surf_usm_h%rad_sw_dif = 0.0_wp 
-          surf_usm_h%rad_sw_ref = 0.0_wp 
-          surf_usm_h%rad_sw_res = 0.0_wp 
-          surf_usm_h%rad_lw_in  = 0.0_wp 
-          surf_usm_h%rad_lw_out = 0.0_wp 
-          surf_usm_h%rad_lw_dif = 0.0_wp 
-          surf_usm_h%rad_lw_ref = 0.0_wp 
-          surf_usm_h%rad_lw_res = 0.0_wp 
+          surf_usm_h%rad_sw_in  = 0.0_wp
+          surf_usm_h%rad_sw_out = 0.0_wp
+          surf_usm_h%rad_sw_dir = 0.0_wp
+          surf_usm_h%rad_sw_dif = 0.0_wp
+          surf_usm_h%rad_sw_ref = 0.0_wp
+          surf_usm_h%rad_sw_res = 0.0_wp
+          surf_usm_h%rad_lw_in  = 0.0_wp
+          surf_usm_h%rad_lw_out = 0.0_wp
+          surf_usm_h%rad_lw_dif = 0.0_wp
+          surf_usm_h%rad_lw_ref = 0.0_wp
+          surf_usm_h%rad_lw_res = 0.0_wp
        ENDIF
        DO  l = 0, 3
@@ -1729 +1732 @@
              ALLOCATE( surf_lsm_v(l)%rad_lw_res(1:surf_lsm_v(l)%ns) )
 
-             surf_lsm_v(l)%rad_sw_in  = 0.0_wp 
+             surf_lsm_v(l)%rad_sw_in  = 0.0_wp
              surf_lsm_v(l)%rad_sw_out = 0.0_wp
              surf_lsm_v(l)%rad_sw_dir = 0.0_wp
@@ -1736 +1739 @@
              surf_lsm_v(l)%rad_sw_res = 0.0_wp
 
-             surf_lsm_v(l)%rad_lw_in  = 0.0_wp 
-             surf_lsm_v(l)%rad_lw_out = 0.0_wp 
-             surf_lsm_v(l)%rad_lw_dif = 0.0_wp 
-             surf_lsm_v(l)%rad_lw_ref = 0.0_wp 
-             surf_lsm_v(l)%rad_lw_res = 0.0_wp 
+             surf_lsm_v(l)%rad_lw_in  = 0.0_wp
+             surf_lsm_v(l)%rad_lw_out = 0.0_wp
+             surf_lsm_v(l)%rad_lw_dif = 0.0_wp
+             surf_lsm_v(l)%rad_lw_ref = 0.0_wp
+             surf_lsm_v(l)%rad_lw_res = 0.0_wp
           ENDIF
           IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_sw_in )  .AND.              &
@@ -1755 +1758 @@
              ALLOCATE( surf_usm_v(l)%rad_lw_ref(1:surf_usm_v(l)%ns) )
              ALLOCATE( surf_usm_v(l)%rad_lw_res(1:surf_usm_v(l)%ns) )
-             surf_usm_v(l)%rad_sw_in  = 0.0_wp 
+             surf_usm_v(l)%rad_sw_in  = 0.0_wp
              surf_usm_v(l)%rad_sw_out = 0.0_wp
              surf_usm_v(l)%rad_sw_dir = 0.0_wp
@@ -1761 +1764 @@
              surf_usm_v(l)%rad_sw_ref = 0.0_wp
              surf_usm_v(l)%rad_sw_res = 0.0_wp
-             surf_usm_v(l)%rad_lw_in  = 0.0_wp 
-             surf_usm_v(l)%rad_lw_out = 0.0_wp 
-             surf_usm_v(l)%rad_lw_dif = 0.0_wp 
-             surf_usm_v(l)%rad_lw_ref = 0.0_wp 
-             surf_usm_v(l)%rad_lw_res = 0.0_wp 
+             surf_usm_v(l)%rad_lw_in  = 0.0_wp
+             surf_usm_v(l)%rad_lw_out = 0.0_wp
+             surf_usm_v(l)%rad_lw_dif = 0.0_wp
+             surf_usm_v(l)%rad_lw_ref = 0.0_wp
+             surf_usm_v(l)%rad_lw_res = 0.0_wp
           ENDIF
        ENDDO
@@ -1855 +1858 @@
           ENDDO
 !
-!--       Level 2 initialization of broadband albedo via given albedo_type. 
-!--       Only if albedo_type is non-zero. In case of urban surface and 
+!--       Level 2 initialization of broadband albedo via given albedo_type.
+!--       Only if albedo_type is non-zero. In case of urban surface and
 !--       input data is read from ASCII file, albedo_type will be zero, so that
-!--       albedo won't be overwritten. 
+!--       albedo won't be overwritten.
           DO  m = 1, surf_lsm_h%ns
              IF ( surf_lsm_h%albedo_type(m,ind_veg_wall) /= 0 )                &
@@ -1931 +1934 @@
                    surf_usm_h%albedo(m,ind_wat_win)   = albedo_pars_f%pars_xy(0,j,i)
                 ENDIF
-             ENDDO 
-!
-!--          Vertical surfaces            
+             ENDDO
+!
+!--          Vertical surfaces
              DO  l = 0, 3
 
@@ -1961 +1964 @@
              ENDDO
 
-          ENDIF  
+          ENDIF
 !
 !--       Read explicit albedo values from building surface pars. If present,
@@ -2051 +2054 @@
 !--       Allocate albedos for short/longwave radiation, horizontal surfaces
 !--       for wall/green/window (USM) or vegetation/pavement/water surfaces
-!--       (LSM). 
+!--       (LSM).
           ALLOCATE ( surf_lsm_h%aldif(1:surf_lsm_h%ns,0:2)       )
           ALLOCATE ( surf_lsm_h%aldir(1:surf_lsm_h%ns,0:2)       )
@@ -2071 +2074 @@
 
 !
-!--       Allocate broadband albedo (temporary for the current radiation 
+!--       Allocate broadband albedo (temporary for the current radiation
 !--       implementations)
           IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )                            &
@@ -2079 +2082 @@
 
 !
-!--       Allocate albedos for short/longwave radiation, vertical surfaces 
+!--       Allocate albedos for short/longwave radiation, vertical surfaces
           DO  l = 0, 3
 
@@ -2111 +2114 @@
           ENDDO
 !
-!--       Level 1 initialization of spectral albedos via namelist 
-!--       paramters. Please note, this case all surface tiles are initialized 
-!--       the same. 
+!--       Level 1 initialization of spectral albedos via namelist
+!--       paramters. Please note, this case all surface tiles are initialized
+!--       the same.
           IF ( surf_lsm_h%ns > 0 )  THEN
              surf_lsm_h%aldif  = albedo_lw_dif
@@ -2162 +2165 @@
 
 !
-!--       Level 2 initialization of spectral albedos via albedo_type. 
-!--       Please note, for natural- and urban-type surfaces, a tile approach 
-!--       is applied so that the resulting albedo is calculated via the weighted 
-!--       average of respective surface fractions. 
+!--       Level 2 initialization of spectral albedos via albedo_type.
+!--       Please note, for natural- and urban-type surfaces, a tile approach
+!--       is applied so that the resulting albedo is calculated via the weighted
+!--       average of respective surface fractions.
           DO  m = 1, surf_lsm_h%ns
 !
@@ -2186 +2189 @@
           ENDDO
 !
-!--       For urban surface only if albedo has not been already initialized 
+!--       For urban surface only if albedo has not been already initialized
 !--       in the urban-surface model via the ASCII file.
           IF ( .NOT. surf_usm_h%albedo_from_ascii )  THEN
@@ -2231 +2234 @@
              ENDDO
 !
-!--          For urban surface only if albedo has not been already initialized 
+!--          For urban surface only if albedo has not been already initialized
 !--          in the urban-surface model via the ASCII file.
              IF ( .NOT. surf_usm_v(l)%albedo_from_ascii )  THEN
@@ -2269 +2272 @@
                    IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )   &
                       surf_lsm_h%albedo(m,ind_type) =                          &
-                                             albedo_pars_f%pars_xy(0,j,i)      
+                                             albedo_pars_f%pars_xy(0,j,i)
                    IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )   &
                       surf_lsm_h%aldir(m,ind_type) =                           &
-                                             albedo_pars_f%pars_xy(1,j,i)      
+                                             albedo_pars_f%pars_xy(1,j,i)
                    IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )   &
                       surf_lsm_h%aldif(m,ind_type) =                           &
-                                             albedo_pars_f%pars_xy(1,j,i)      
+                                             albedo_pars_f%pars_xy(1,j,i)
                    IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )   &
                       surf_lsm_h%asdir(m,ind_type) =                           &
-                                             albedo_pars_f%pars_xy(2,j,i)      
+                                             albedo_pars_f%pars_xy(2,j,i)
                    IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )   &
                       surf_lsm_h%asdif(m,ind_type) =                           &
@@ -2285 +2288 @@
              ENDDO
 !
-!--          For urban surface only if albedo has not been already initialized 
+!--          For urban surface only if albedo has not been already initialized
 !--          in the urban-surface model via the ASCII file.
              IF ( .NOT. surf_usm_h%albedo_from_ascii )  THEN
@@ -2378 +2381 @@
                 ENDDO
 !
-!--             For urban surface only if albedo has not been already initialized 
+!--             For urban surface only if albedo has not been already initialized
 !--             in the urban-surface model via the ASCII file.
                 IF ( .NOT. surf_usm_v(l)%albedo_from_ascii )  THEN
@@ -2411 +2414 @@
                                          albedo_pars_f%pars_xy(2,j+joff,i+ioff)
                       ENDIF
-!                     
+!
 !--                   Spectral albedos especially for building green surfaces
                       IF ( albedo_pars_f%pars_xy(3,j+joff,i+ioff) /=           &
@@ -2427 +2430 @@
                                          albedo_pars_f%pars_xy(4,j+joff,i+ioff)
                       ENDIF
-!                     
+!
 !--                   Spectral albedos especially for building window surfaces
                       IF ( albedo_pars_f%pars_xy(5,j+joff,i+ioff) /=           &
@@ -2640 +2643 @@
 
 !
-!--       Calculate initial values of current (cosine of) the zenith angle and 
+!--       Calculate initial values of current (cosine of) the zenith angle and
 !--       whether the sun is up
           CALL get_date_time( time_since_reference_point, &
@@ -2785 +2788 @@
 
 !
-!--       Allocate 1-element array for surface temperature 
+!--       Allocate 1-element array for surface temperature
 !--       (RRTMG anticipates an array as passed argument).
           ALLOCATE ( rrtm_tsfc(1) )
 !
-!--       Allocate surface emissivity. 
+!--       Allocate surface emissivity.
 !--       Values will be given directly before calling rrtm_lw.
           ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
@@ -2801 +2804 @@
                             '&Please provide <jobname>_lsw file in the INPUT directory.'
              CALL message( 'radiation_init', 'PA0583', 1, 2, 0, 6, 0 )
-          ENDIF         
+          ENDIF
           INQUIRE( FILE='rrtmg_sw.nc', EXIST=sw_exists )
           IF ( .NOT. sw_exists )  THEN
@@ -2808 +2811 @@
                             '&Please provide <jobname>_rsw file in the INPUT directory.'
              CALL message( 'radiation_init', 'PA0584', 1, 2, 0, 6, 0 )
-          ENDIF          
-          
+          ENDIF
+
           IF ( lw_radiation )  CALL rrtmg_lw_ini ( c_p )
           IF ( sw_radiation )  CALL rrtmg_sw_ini ( c_p )
-          
+
 !
 !--       Set input files for RRTMG
-          INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists) 
+          INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists)
           IF ( .NOT. snd_exists )  THEN
              rrtm_input_file = "rrtmg_lw.nc"
@@ -2823 +2826 @@
 !--       Read vertical layers for RRTMG from sounding data
 !--       The routine provides nzt_rad, hyp_snd(1:nzt_rad),
-!--       t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1), 
+!--       t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
 !--       rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
           CALL read_sounding_data
@@ -2834 +2837 @@
        ENDIF
 !
-!--    Initializaion actions exclusively required for external 
+!--    Initializaion actions exclusively required for external
 !--    radiation forcing
        IF ( radiation_scheme == 'external' )  THEN
 !
-!--       Open the radiation input file. Note, for child domain, a dynamic 
-!--       input file is often not provided. In order to do not need to 
+!--       Open the radiation input file. Note, for child domain, a dynamic
+!--       input file is often not provided. In order to do not need to
 !--       duplicate the dynamic input file just for the radiation input, take
-!--       it from the dynamic file for the parent if not available for the 
-!--       child domain(s). In this case this is possible because radiation 
+!--       it from the dynamic file for the parent if not available for the
+!--       child domain(s). In this case this is possible because radiation
 !--       input should be the same for each model.
           INQUIRE( FILE = TRIM( input_file_dynamic ),                          &
                    EXIST = radiation_input_root_domain  )
-                   
+
           IF ( .NOT. input_pids_dynamic  .AND.                                 &
                .NOT. radiation_input_root_domain )  THEN
@@ -2859 +2862 @@
 !
 !--       Open dynamic input file for child domain if available, else, open
-!--       dynamic input file for the root domain. 
+!--       dynamic input file for the root domain.
           IF ( input_pids_dynamic )  THEN
              CALL open_read_file( TRIM( input_file_dynamic ) //                &
@@ -2868 +2871 @@
                                   pids_id )
           ENDIF
-                               
+
           CALL inquire_num_variables( pids_id, num_var_pids )
-!         
+!
 !--       Allocate memory to store variable names and read them
           ALLOCATE( vars_pids(1:num_var_pids) )
           CALL inquire_variable_names( pids_id, vars_pids )
-!         
+!
 !--       Input time dimension.
           IF ( check_existence( vars_pids, 'time_rad' ) )  THEN
              CALL get_dimension_length( pids_id, ntime, 'time_rad' )
-          
+
              ALLOCATE( time_rad_f%var1d(0:ntime-1) )
-!                                                                                 
-!--          Read variable                    
+!
+!--          Read variable
              CALL get_variable( pids_id, 'time_rad', time_rad_f%var1d )
-                                
+
              time_rad_f%from_file = .TRUE.
-          ENDIF           
-!         
+          ENDIF
+!
 !--       Input shortwave downwelling.
           IF ( check_existence( vars_pids, 'rad_sw_in' ) )  THEN
-!         
+!
 !--          Get _FillValue attribute
              CALL get_attribute( pids_id, char_fill, rad_sw_in_f%fill,         &
                                  .FALSE., 'rad_sw_in' )
-!         
+!
 !--          Get level-of-detail
              CALL get_attribute( pids_id, char_lod, rad_sw_in_f%lod,           &
@@ -2905 +2908 @@
 !
 !--          Level-of-detail 2 - radiation depends on time_rad, y, x
-             ELSEIF ( rad_sw_in_f%lod == 2 )  THEN  
+             ELSEIF ( rad_sw_in_f%lod == 2 )  THEN
                 ALLOCATE( rad_sw_in_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )
-                
+
                 CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var3d,    &
                                    nxl, nxr, nys, nyn, 0, ntime-1 )
-                                   
+
                 rad_sw_in_f%from_file = .TRUE.
-             ELSE 
+             ELSE
                 message_string = '"rad_sw_in" has no valid lod attribute'
                 CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
              ENDIF
           ENDIF
-!         
+!
 !--       Input longwave downwelling.
           IF ( check_existence( vars_pids, 'rad_lw_in' ) )  THEN
-!         
+!
 !--          Get _FillValue attribute
              CALL get_attribute( pids_id, char_fill, rad_lw_in_f%fill,         &
                                  .FALSE., 'rad_lw_in' )
-!         
+!
 !--          Get level-of-detail
              CALL get_attribute( pids_id, char_lod, rad_lw_in_f%lod,           &
@@ -2936 +2939 @@
 !
 !--          Level-of-detail 2 - radiation depends on time_rad, y, x
-             ELSEIF ( rad_lw_in_f%lod == 2 )  THEN  
+             ELSEIF ( rad_lw_in_f%lod == 2 )  THEN
                 ALLOCATE( rad_lw_in_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )
-                
+
                 CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var3d,    &
                                    nxl, nxr, nys, nyn, 0, ntime-1 )
-                                   
+
                 rad_lw_in_f%from_file = .TRUE.
-             ELSE 
+             ELSE
                 message_string = '"rad_lw_in" has no valid lod attribute'
                 CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
              ENDIF
           ENDIF
-!         
+!
 !--       Input shortwave downwelling, diffuse part.
           IF ( check_existence( vars_pids, 'rad_sw_in_dif' ) )  THEN
-!         
+!
 !--          Read _FillValue attribute
              CALL get_attribute( pids_id, char_fill, rad_sw_in_dif_f%fill,     &
                                  .FALSE., 'rad_sw_in_dif' )
-!         
+!
 !--          Get level-of-detail
              CALL get_attribute( pids_id, char_lod, rad_sw_in_dif_f%lod,       &
@@ -2968 +2971 @@
 !
 !--          Level-of-detail 2 - radiation depends on time_rad, y, x
-             ELSEIF ( rad_sw_in_dif_f%lod == 2 )  THEN  
+             ELSEIF ( rad_sw_in_dif_f%lod == 2 )  THEN
                 ALLOCATE( rad_sw_in_dif_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )
-                
+
                 CALL get_variable( pids_id, 'rad_sw_in_dif',                   &
                                    rad_sw_in_dif_f%var3d,                      &
                                    nxl, nxr, nys, nyn, 0, ntime-1 )
-                                   
+
                 rad_sw_in_dif_f%from_file = .TRUE.
-             ELSE 
+             ELSE
                 message_string = '"rad_sw_in_dif" has no valid lod attribute'
                 CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
              ENDIF
           ENDIF
-!         
+!
 !--       Finally, close the input file and deallocate temporary arrays
           DEALLOCATE( vars_pids )
-          
+
           CALL close_input_file( pids_id )
 #endif
@@ -2995 +2998 @@
              CALL message( 'radiation_init', 'PA0195', 1, 2, 0, 6, 0 )
           ENDIF
-          
+
           IF ( .NOT. time_rad_f%from_file )  THEN
              message_string = 'In case of external radiation forcing ' //      &
@@ -3001 +3004 @@
              CALL message( 'radiation_init', 'PA0196', 1, 2, 0, 6, 0 )
           ENDIF
-          
+
           CALL get_date_time( 0.0_wp, second_of_day=second_of_day )
 
@@ -3018 +3021 @@
              ENDIF
           ENDIF
-          
+
           IF ( ALLOCATED( rad_lw_in_f%var1d ) )  THEN
              IF ( ANY( rad_lw_in_f%var1d == rad_lw_in_f%fill ) )  THEN
@@ -3026 +3029 @@
              ENDIF
           ENDIF
-          
+
           IF ( ALLOCATED( rad_sw_in_dif_f%var1d ) )  THEN
              IF ( ANY( rad_sw_in_dif_f%var1d == rad_sw_in_dif_f%fill ) )  THEN
@@ -3034 +3037 @@
              ENDIF
           ENDIF
-          
+
           IF ( ALLOCATED( rad_sw_in_f%var3d ) )  THEN
              IF ( ANY( rad_sw_in_f%var3d == rad_sw_in_f%fill ) )  THEN
@@ -3042 +3045 @@
              ENDIF
           ENDIF
-          
+
           IF ( ALLOCATED( rad_lw_in_f%var3d ) )  THEN
              IF ( ANY( rad_lw_in_f%var3d == rad_lw_in_f%fill ) )  THEN
@@ -3050 +3053 @@
              ENDIF
           ENDIF
-          
+
           IF ( ALLOCATED( rad_sw_in_dif_f%var3d ) )  THEN
              IF ( ANY( rad_sw_in_dif_f%var3d == rad_sw_in_dif_f%fill ) )  THEN
@@ -3068 +3071 @@
           ENDIF
 !
-!--       All radiation input should have the same level of detail. The sum 
+!--       All radiation input should have the same level of detail. The sum
 !--       of lods divided by the number of available radiation arrays must be
-!--       1 (if all are lod = 1) or 2 (if all are lod = 2). 
+!--       1 (if all are lod = 1) or 2 (if all are lod = 2).
           IF ( REAL( MERGE( rad_lw_in_f%lod, 0, rad_lw_in_f%from_file ) +       &
                      MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) +       &
@@ -3109 +3112 @@
           CASE ( 'constant' )
              CALL radiation_constant
-             
+
           CASE ( 'external' )
 !
@@ -3128 +3131 @@
 
 !
-!--    If required, read or calculate and write out the SVF 
+!--    If required, read or calculate and write out the SVF
        IF ( radiation_interactions .AND. read_svf)  THEN
 !
@@ -3147 +3150 @@
 
 !
-!--    Adjust radiative fluxes. In case of urban and land surfaces, also 
+!--    Adjust radiative fluxes. In case of urban and land surfaces, also
 !--    call an initial interaction.
        IF ( radiation_interactions )  THEN
@@ -3171 +3174 @@
        INTEGER(iwp) ::  l   !< running index for surface orientation
        INTEGER(iwp) ::  t   !< index of current timestep
-       INTEGER(iwp) ::  tm  !< index of previous timestep 
-       
+       INTEGER(iwp) ::  tm  !< index of previous timestep
+
        LOGICAL      ::  horizontal !< flag indicating treatment of horinzontal surfaces
-       
-       REAL(wp) ::  fac_dt               !< interpolation factor  
+
+       REAL(wp) ::  fac_dt               !< interpolation factor
        REAL(wp) ::  second_of_day_init   !< second of the day at model start
 
-       TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine   
+       TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine
 
 !
@@ -3187 +3190 @@
        CALL calc_zenith( day_of_year, second_of_day )
 !
-!--    Interpolate external radiation on current timestep 
+!--    Interpolate external radiation on current timestep
        IF ( time_since_reference_point  <= 0.0_wp )  THEN
           t      = 0
@@ -3198 +3201 @@
              t = t + 1
           ENDDO
-          
+
           tm = MAX( t-1, 0 )
-          
+
           fac_dt = ( time_since_reference_point                                &
                    - time_rad_f%var1d(tm) + dt_3d )                            &
@@ -3207 +3210 @@
        ENDIF
 !
-!--    Call clear-sky calculation for each surface orientation. 
+!--    Call clear-sky calculation for each surface orientation.
 !--    First, horizontal surfaces
        horizontal = .TRUE.
@@ -3223 +3226 @@
           CALL radiation_external_surf
        ENDDO
-       
+
        CONTAINS
 
@@ -3229 +3232 @@
 
              USE control_parameters
-          
+
              IMPLICIT NONE
 
-             INTEGER(iwp) ::  i    !< grid index along x-dimension   
-             INTEGER(iwp) ::  j    !< grid index along y-dimension  
-             INTEGER(iwp) ::  k    !< grid index along z-dimension    
+             INTEGER(iwp) ::  i    !< grid index along x-dimension
+             INTEGER(iwp) ::  j    !< grid index along y-dimension
+             INTEGER(iwp) ::  k    !< grid index along z-dimension
              INTEGER(iwp) ::  m    !< running index for surface elements
-             
-             REAL(wp) ::  lw_in     !< downwelling longwave radiation, interpolated value     
+
+             REAL(wp) ::  lw_in     !< downwelling longwave radiation, interpolated value
              REAL(wp) ::  sw_in     !< downwelling shortwave radiation, interpolated value
-             REAL(wp) ::  sw_in_dif !< downwelling diffuse shortwave radiation, interpolated value   
+             REAL(wp) ::  sw_in_dif !< downwelling diffuse shortwave radiation, interpolated value
 
              IF ( surf%ns < 1 )  RETURN
@@ -3245 +3248 @@
 !--          level-of-detail = 1. Note, here it must be distinguished between
 !--          averaged radiation and non-averaged radiation for the upwelling
-!--          fluxes. 
+!--          fluxes.
              IF ( rad_sw_in_f%lod == 1 )  THEN
-             
+
                 sw_in = ( 1.0_wp - fac_dt ) * rad_sw_in_f%var1d(tm)            &
                                    + fac_dt * rad_sw_in_f%var1d(t)
-                                         
+
                 lw_in = ( 1.0_wp - fac_dt ) * rad_lw_in_f%var1d(tm)            &
                                    + fac_dt * rad_lw_in_f%var1d(t)
 !
-!--             Limit shortwave incoming radiation to positive values, in order 
+!--             Limit shortwave incoming radiation to positive values, in order
 !--             to overcome possible observation errors.
                 sw_in = MAX( 0.0_wp, sw_in )
                 sw_in = MERGE( sw_in, 0.0_wp, sun_up )
-                          
-                surf%rad_sw_in = sw_in                                         
+
+                surf%rad_sw_in = sw_in
                 surf%rad_lw_in = lw_in
-             
+
                 IF ( average_radiation )  THEN
                    surf%rad_sw_out = albedo_urb * surf%rad_sw_in
-                   
+
                    surf%rad_lw_out = emissivity_urb * sigma_sb * t_rad_urb**4  &
                                   + ( 1.0_wp - emissivity_urb ) * surf%rad_lw_in
-                   
+
                    surf%rad_net = surf%rad_sw_in - surf%rad_sw_out             &
                                 + surf%rad_lw_in - surf%rad_lw_out
-                   
+
                    surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb          &
                                                      * sigma_sb                &
@@ -3284 +3287 @@
                                              surf%albedo(m,ind_wat_win) )      &
                                            * surf%rad_sw_in(m)
-                   
+
                       surf%rad_lw_out(m) = ( surf%frac(m,ind_veg_wall)  *      &
                                              surf%emissivity(m,ind_veg_wall)   &
@@ -3294 +3297 @@
                                            * sigma_sb                          &
                                            * ( surf%pt_surface(m) * exner(k) )**4
-                   
+
                       surf%rad_lw_out_change_0(m) =                            &
                                          ( surf%frac(m,ind_veg_wall)  *        &
@@ -3305 +3308 @@
                                          * ( surf%pt_surface(m) * exner(k) )**3
                    ENDDO
-                
+
                 ENDIF
 !
-!--             If diffuse shortwave radiation is available, store it on 
+!--             If diffuse shortwave radiation is available, store it on
 !--             the respective files.
-                IF ( rad_sw_in_dif_f%from_file )  THEN 
+                IF ( rad_sw_in_dif_f%from_file )  THEN
                    sw_in_dif= ( 1.0_wp - fac_dt ) * rad_sw_in_dif_f%var1d(tm)  &
                                          + fac_dt * rad_sw_in_dif_f%var1d(t)
-                                         
+
                    IF ( ALLOCATED( rad_sw_in_diff ) )  rad_sw_in_diff = sw_in_dif
                    IF ( ALLOCATED( rad_sw_in_dir  ) )  rad_sw_in_dir  = sw_in  &
                                                                     - sw_in_dif
-!              
-!--                Diffuse longwave radiation equals the total downwelling 
+!
+!--                Diffuse longwave radiation equals the total downwelling
 !--                longwave radiation
                    IF ( ALLOCATED( rad_lw_in_diff ) )  rad_lw_in_diff = lw_in
@@ -3330 +3333 @@
                    j = surf%j(m)
                    k = surf%k(m)
-                   
+
                    surf%rad_sw_in(m) = ( 1.0_wp - fac_dt )                     &
                                             * rad_sw_in_f%var3d(tm,j,i)        &
-                                   + fac_dt * rad_sw_in_f%var3d(t,j,i)                                   
-!
-!--                Limit shortwave incoming radiation to positive values, in 
+                                   + fac_dt * rad_sw_in_f%var3d(t,j,i)
+!
+!--                Limit shortwave incoming radiation to positive values, in
 !--                order to overcome possible observation errors.
                    surf%rad_sw_in(m) = MAX( 0.0_wp, surf%rad_sw_in(m) )
                    surf%rad_sw_in(m) = MERGE( surf%rad_sw_in(m), 0.0_wp, sun_up )
-                                         
+
                    surf%rad_lw_in(m) = ( 1.0_wp - fac_dt )                     &
                                             * rad_lw_in_f%var3d(tm,j,i)        &
-                                   + fac_dt * rad_lw_in_f%var3d(t,j,i) 
-!
-!--                Weighted average according to surface fraction. 
+                                   + fac_dt * rad_lw_in_f%var3d(t,j,i)
+!
+!--                Weighted average according to surface fraction.
                    surf%rad_sw_out(m) = ( surf%frac(m,ind_veg_wall)  *         &
                                           surf%albedo(m,ind_veg_wall)          &
@@ -3376 +3379 @@
                                    + surf%rad_lw_in(m) - surf%rad_lw_out(m)
 !
-!--                If diffuse shortwave radiation is available, store it on 
-!--                the respective files. 
+!--                If diffuse shortwave radiation is available, store it on
+!--                the respective files.
                    IF ( rad_sw_in_dif_f%from_file )  THEN
                       IF ( ALLOCATED( rad_sw_in_diff ) )                       &
@@ -3384 +3387 @@
                                      + fac_dt * rad_sw_in_dif_f%var3d(t,j,i)
 !
-!--                   dir = sw_in - sw_in_dif. 
+!--                   dir = sw_in - sw_in_dif.
                       IF ( ALLOCATED( rad_sw_in_dir  ) )                       &
                          rad_sw_in_dir(j,i)  = surf%rad_sw_in(m) -             &
                                                rad_sw_in_diff(j,i)
-!                 
-!--                   Diffuse longwave radiation equals the total downwelling 
+!
+!--                   Diffuse longwave radiation equals the total downwelling
 !--                   longwave radiation
                       IF ( ALLOCATED( rad_lw_in_diff ) )                       &
@@ -3399 +3402 @@
              ENDIF
 !
-!--          Store radiation also on 2D arrays, which are still used for 
+!--          Store radiation also on 2D arrays, which are still used for
 !--          direct-diffuse splitting. Note, this is only required
 !--          for horizontal surfaces, which covers all x,y position.
@@ -3406 +3409 @@
                    i = surf%i(m)
                    j = surf%j(m)
-                   
+
                    rad_sw_in(0,j,i)  = surf%rad_sw_in(m)
                    rad_lw_in(0,j,i)  = surf%rad_lw_in(m)
@@ -3413 +3416 @@
                 ENDDO
              ENDIF
- 
+
           END SUBROUTINE radiation_external_surf
 
-    END SUBROUTINE radiation_external    
-    
+    END SUBROUTINE radiation_external
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -3431 +3434 @@
        LOGICAL      ::  horizontal !< flag indicating treatment of horinzontal surfaces
 
-       REAL(wp)     ::  pt1       !< potential temperature at first grid level or mean value at urban layer top 
-       REAL(wp)     ::  pt1_l     !< potential temperature at first grid level or mean value at urban layer top at local subdomain 
-       REAL(wp)     ::  ql1       !< liquid water mixing ratio at first grid level or mean value at urban layer top 
-       REAL(wp)     ::  ql1_l     !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain 
-
-       TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine   
+       REAL(wp)     ::  pt1       !< potential temperature at first grid level or mean value at urban layer top
+       REAL(wp)     ::  pt1_l     !< potential temperature at first grid level or mean value at urban layer top at local subdomain
+       REAL(wp)     ::  ql1       !< liquid water mixing ratio at first grid level or mean value at urban layer top
+       REAL(wp)     ::  ql1_l     !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain
+
+       TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine
 
 !
@@ -3452 +3455 @@
 !--    Calculate value of the Exner function at model surface
 !
-!--    In case averaged radiation is used, calculate mean temperature and 
+!--    In case averaged radiation is used, calculate mean temperature and
 !--    liquid water mixing ratio at the urban-layer top.
        IF ( average_radiation ) THEN
@@ -3461 +3464 @@
           IF ( bulk_cloud_model  .OR.  cloud_droplets  )  ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) )
 
-#if defined( __parallel )      
+#if defined( __parallel )
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
           CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
@@ -3477 +3480 @@
           ENDIF
 #else
-          pt1 = pt1_l 
+          pt1 = pt1_l
           IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1 = ql1_l
 #endif
@@ -3487 +3490 @@
        ENDIF
 !
-!--    Call clear-sky calculation for each surface orientation. 
+!--    Call clear-sky calculation for each surface orientation.
 !--    First, horizontal surfaces
        horizontal = .TRUE.
@@ -3518 +3521 @@
 
 !
-!--          Calculate radiation fluxes and net radiation (rad_net) assuming 
-!--          homogeneous urban radiation conditions. 
-             IF ( average_radiation ) THEN       
+!--          Calculate radiation fluxes and net radiation (rad_net) assuming
+!--          homogeneous urban radiation conditions.
+             IF ( average_radiation ) THEN
 
                 k = nz_urban_t
@@ -3526 +3529 @@
                 surf%rad_sw_in  = solar_constant * sky_trans * cos_zenith
                 surf%rad_sw_out = albedo_urb * surf%rad_sw_in
-                
+
                 surf%rad_lw_in  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k+1))**4
 
@@ -3539 +3542 @@
 
 !
-!--          Calculate radiation fluxes and net radiation (rad_net) for each surface 
+!--          Calculate radiation fluxes and net radiation (rad_net) for each surface
 !--          element.
              ELSE
@@ -3551 +3554 @@
 
 !
-!--                Weighted average according to surface fraction. 
-!--                ATTENTION: when radiation interactions are switched on the 
-!--                calculated fluxes below are not actually used as they are 
+!--                Weighted average according to surface fraction.
+!--                ATTENTION: when radiation interactions are switched on the
+!--                calculated fluxes below are not actually used as they are
 !--                overwritten in radiation_interaction.
                    surf%rad_sw_out(m) = ( surf%frac(m,ind_veg_wall)  *         &
@@ -3611 +3614 @@
                 ENDDO
              ENDIF
- 
+
           END SUBROUTINE radiation_clearsky_surf
 
@@ -3628 +3631 @@
 
        INTEGER(iwp) ::  l         !< running index for surface orientation
-       
+
        LOGICAL      ::  horizontal !< flag indicating treatment of horinzontal surfaces
 
-       REAL(wp)     ::  pt1       !< potential temperature at first grid level or mean value at urban layer top 
-       REAL(wp)     ::  pt1_l     !< potential temperature at first grid level or mean value at urban layer top at local subdomain 
-       REAL(wp)     ::  ql1       !< liquid water mixing ratio at first grid level or mean value at urban layer top 
-       REAL(wp)     ::  ql1_l     !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain 
+       REAL(wp)     ::  pt1       !< potential temperature at first grid level or mean value at urban layer top
+       REAL(wp)     ::  pt1_l     !< potential temperature at first grid level or mean value at urban layer top at local subdomain
+       REAL(wp)     ::  ql1       !< liquid water mixing ratio at first grid level or mean value at urban layer top
+       REAL(wp)     ::  ql1_l     !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain
 
        TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine
 
 !
-!--    In case averaged radiation is used, calculate mean temperature and 
+!--    In case averaged radiation is used, calculate mean temperature and
 !--    liquid water mixing ratio at the urban-layer top.
-       IF ( average_radiation ) THEN   
+       IF ( average_radiation ) THEN
           pt1   = 0.0_wp
           IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1   = 0.0_wp
@@ -3648 +3651 @@
           IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) )
 
-#if defined( __parallel )      
+#if defined( __parallel )
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
           CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
@@ -3729 +3732 @@
              ELSE
 !
-!--             Determine index offset between surface element and adjacent 
+!--             Determine index offset between surface element and adjacent
 !--             atmospheric grid point
                 ioff = surf%ioff
@@ -3753 +3756 @@
 
 !
-!--                Weighted average according to surface fraction. 
+!--                Weighted average according to surface fraction.
                    surf%rad_lw_out(m) = ( surf%frac(m,ind_veg_wall)  *         &
                                           surf%emissivity(m,ind_veg_wall)      &
@@ -3802 +3805 @@
 
           END SUBROUTINE radiation_constant_surf
-          
+
 
     END SUBROUTINE radiation_constant
@@ -3815 +3818 @@
 
        IMPLICIT NONE
- 
+
        INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
-    
-
-       
+
+
+
 !
 !--    Write radiation model header
@@ -3834 +3837 @@
        ELSEIF ( radiation_scheme == "external" )  THEN
           WRITE( io, 14 )
-       ENDIF 
+       ENDIF
 
        IF ( albedo_type_f%from_file  .OR.  vegetation_type_f%from_file  .OR.   &
@@ -3840 +3843 @@
             building_type_f%from_file )  THEN
              WRITE( io, 13 )
-       ELSE  
+       ELSE
           IF ( albedo_type == 0 )  THEN
              WRITE( io, 7 ) albedo
@@ -3850 +3853 @@
           WRITE( io, 9 )
        ENDIF
-       
+
        WRITE( io, 12 ) dt_radiation
- 
+
 
  3 FORMAT (//' Radiation model information:'/                                  &
@@ -3873 +3876 @@
 
     END SUBROUTINE radiation_header
-   
+
 
 !------------------------------------------------------------------------------!
@@ -3885 +3888 @@
        IMPLICIT NONE
 
-       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file 
-       
+       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
+
        NAMELIST /radiation_par/   albedo, albedo_lw_dif, albedo_lw_dir,         &
                                   albedo_sw_dif, albedo_sw_dir, albedo_type,    &
@@ -3903 +3906 @@
                                   unscheduled_radiation_calls
 
-   
+
        NAMELIST /radiation_parameters/ albedo, albedo_lw_dif, albedo_lw_dir,    &
                                   albedo_sw_dif, albedo_sw_dir, albedo_type,    &
@@ -3919 +3922 @@
                                   svfnorm_report_thresh, sw_radiation,          &
                                   unscheduled_radiation_calls
-   
+
        line = ' '
-       
+
 !
 !--    Try to find radiation model namelist
@@ -3979 +3982 @@
 
  14    CONTINUE
-       
+
     END SUBROUTINE radiation_parin
 
@@ -4039 +4042 @@
        zenith(0) = cos_zenith
 !
-!--    Calculate surface albedo. In case average radiation is applied, 
+!--    Calculate surface albedo. In case average radiation is applied,
 !--    this is not required.
 #if defined( __netcdf )
@@ -4071 +4074 @@
        IF ( average_radiation ) THEN
 !
-!--       Determine minimum topography top index. 
+!--       Determine minimum topography top index.
           k_topo_l = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
 #if defined( __parallel )
@@ -4079 +4082 @@
           k_topo = k_topo_l
 #endif
-       
+
           rrtm_asdir(1)  = albedo_urb
           rrtm_asdif(1)  = albedo_urb
@@ -4087 +4090 @@
           rrtm_emis = emissivity_urb
 !
-!--       Calculate mean pt profile. 
+!--       Calculate mean pt profile.
           CALL calc_mean_profile( pt, 4 )
           pt_av = hom(:, 1, 4, 0)
-          
+
           IF ( humidity )  THEN
              CALL calc_mean_profile( q, 41 )
@@ -4104 +4107 @@
              ! average ql is now in hom(:, 1, 54, 0)
              ql_av = hom(:, 1, 54, 0)
-             
+
              DO k = nzb+1, nzt+1
                 rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp       &
@@ -4126 +4129 @@
 
 !
-!--       Avoid temperature/humidity jumps at the top of the PALM domain by 
+!--       Avoid temperature/humidity jumps at the top of the PALM domain by
 !--       linear interpolation from nzt+2 to nzt+7. Jumps are induced by
 !--       discrepancies between the values in the  domain and those above that
@@ -4163 +4166 @@
                 rrtm_cliqwp(0,k) =  ql_av(k) * 1000._wp *                   &
                                     (rrtm_plev(0,k) - rrtm_plev(0,k+1))     &
-                                    * 100._wp / g 
+                                    * 100._wp / g
 
                 IF ( rrtm_cliqwp(0,k) > 0._wp )  THEN
@@ -4189 +4192 @@
 !--       Set surface temperature
           rrtm_tsfc = t_rad_urb
-          
-          IF ( lw_radiation )  THEN  
+
+          IF ( lw_radiation )  THEN
 !
 !--          Due to technical reasons, copy optical depth to dummy arguments
 !--          which are allocated on the exact size as the rrtmg_lw is called.
 !--          As one dimesion is allocated with zero size, compiler complains
-!--          that rank of the array does not match that of the 
-!--          assumed-shaped arguments in the RRTMG library. In order to 
-!--          avoid this, write to dummy arguments and give pass the entire 
-!--          dummy array. Seems to be the only existing work-around.  
+!--          that rank of the array does not match that of the
+!--          assumed-shaped arguments in the RRTMG library. In order to
+!--          avoid this, write to dummy arguments and give pass the entire
+!--          dummy array. Seems to be the only existing work-around.
              ALLOCATE( rrtm_lw_taucld_dum(1:nbndlw+1,0:0,k_topo+1:nzt_rad+1) )
              ALLOCATE( rrtm_lw_tauaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1) )
@@ -4206 +4209 @@
              rrtm_lw_tauaer_dum =                                              &
                              rrtm_lw_tauaer(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1)
-          
-             CALL rrtmg_lw( 1,                                                 &                                        
+
+             CALL rrtmg_lw( 1,                                                 &
                             nzt_rad-k_topo,                                    &
                             rrtm_icld,                                         &
@@ -4234 +4237 @@
                             rrtm_cicewp(:,k_topo+1:),                 &
                             rrtm_cliqwp(:,k_topo+1:),                 &
-                            rrtm_reice(:,k_topo+1:),                  & 
+                            rrtm_reice(:,k_topo+1:),                  &
                             rrtm_reliq(:,k_topo+1:),                  &
                             rrtm_lw_tauaer_dum,                                &
@@ -4245 +4248 @@
                             rrtm_lwuflx_dt(:,k_topo:),                &
                             rrtm_lwuflxc_dt(:,k_topo:) )
-                            
+
              DEALLOCATE ( rrtm_lw_taucld_dum )
              DEALLOCATE ( rrtm_lw_tauaer_dum )
@@ -4256 +4259 @@
              rad_lw_in_diff(:,:) = rad_lw_in(k_topo,:,:)
 !
-!--          Save heating rates (convert from K/d to K/h). 
+!--          Save heating rates (convert from K/d to K/h).
 !--          Further, even though an aggregated radiation is computed, map
-!--          signle-column profiles on top of any topography, in order to 
+!--          signle-column profiles on top of any topography, in order to
 !--          obtain correct near surface radiation heating/cooling rates.
              DO  i = nxl, nxr
@@ -4274 +4277 @@
           IF ( sw_radiation .AND. sun_up )  THEN
 !
-!--          Due to technical reasons, copy optical depths and other 
-!--          to dummy arguments which are allocated on the exact size as the 
+!--          Due to technical reasons, copy optical depths and other
+!--          to dummy arguments which are allocated on the exact size as the
 !--          rrtmg_sw is called.
 !--          As one dimesion is allocated with zero size, compiler complains
-!--          that rank of the array does not match that of the 
-!--          assumed-shaped arguments in the RRTMG library. In order to 
-!--          avoid this, write to dummy arguments and give pass the entire 
-!--          dummy array. Seems to be the only existing work-around.  
+!--          that rank of the array does not match that of the
+!--          assumed-shaped arguments in the RRTMG library. In order to
+!--          avoid this, write to dummy arguments and give pass the entire
+!--          dummy array. Seems to be the only existing work-around.
              ALLOCATE( rrtm_sw_taucld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
              ALLOCATE( rrtm_sw_ssacld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
@@ -4290 +4293 @@
              ALLOCATE( rrtm_sw_asmaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
              ALLOCATE( rrtm_sw_ecaer_dum(0:0,k_topo+1:nzt_rad+1,1:naerec+1)  )
-     
+
              rrtm_sw_taucld_dum = rrtm_sw_taucld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
              rrtm_sw_ssacld_dum = rrtm_sw_ssacld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
@@ -4309 +4312 @@
                             rrtm_tlev(:,k_topo+1:nzt_rad+2),                   &
                             rrtm_tsfc,                                         &
-                            rrtm_h2ovmr(:,k_topo+1:nzt_rad+1),                 &                               
-                            rrtm_o3vmr(:,k_topo+1:nzt_rad+1),                  &        
+                            rrtm_h2ovmr(:,k_topo+1:nzt_rad+1),                 &
+                            rrtm_o3vmr(:,k_topo+1:nzt_rad+1),                  &
                             rrtm_co2vmr(:,k_topo+1:nzt_rad+1),                 &
                             rrtm_ch4vmr(:,k_topo+1:nzt_rad+1),                 &
                             rrtm_n2ovmr(:,k_topo+1:nzt_rad+1),                 &
                             rrtm_o2vmr(:,k_topo+1:nzt_rad+1),                  &
-                            rrtm_asdir,                                        & 
+                            rrtm_asdir,                                        &
                             rrtm_asdif,                                        &
                             rrtm_aldir,                                        &
@@ -4339 +4342 @@
                             rrtm_sw_asmaer_dum,                                &
                             rrtm_sw_ecaer_dum,                                 &
-                            rrtm_swuflx(:,k_topo:nzt_rad+1),                   &  
-                            rrtm_swdflx(:,k_topo:nzt_rad+1),                   & 
-                            rrtm_swhr(:,k_topo+1:nzt_rad+1),                   & 
-                            rrtm_swuflxc(:,k_topo:nzt_rad+1),                  & 
+                            rrtm_swuflx(:,k_topo:nzt_rad+1),                   &
+                            rrtm_swdflx(:,k_topo:nzt_rad+1),                   &
+                            rrtm_swhr(:,k_topo+1:nzt_rad+1),                   &
+                            rrtm_swuflxc(:,k_topo:nzt_rad+1),                  &
                             rrtm_swdflxc(:,k_topo:nzt_rad+1),                  &
                             rrtm_swhrc(:,k_topo+1:nzt_rad+1),                  &
                             rrtm_dirdflux(:,k_topo:nzt_rad+1),                 &
                             rrtm_difdflux(:,k_topo:nzt_rad+1) )
-                            
+
              DEALLOCATE( rrtm_sw_taucld_dum )
              DEALLOCATE( rrtm_sw_ssacld_dum )
@@ -4356 +4359 @@
              DEALLOCATE( rrtm_sw_asmaer_dum )
              DEALLOCATE( rrtm_sw_ecaer_dum )
-  
+
 !
 !--          Save radiation fluxes for the entire depth of the model domain
@@ -4382 +4385 @@
           ENDIF
 !
-!--    RRTMG is called for each (j,i) grid point separately, starting at the 
+!--    RRTMG is called for each (j,i) grid point separately, starting at the
 !--    highest topography level. Here no RTM is used since average_radiation is false
        ELSE
@@ -4410 +4413 @@
                       rrtm_tlay(0,k) = pt(k,j,i) * exner(k)                     &
                                         + lv_d_cp * ql(k,j,i)
-                      rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q(k,j,i) 
+                      rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q(k,j,i)
                    ENDDO
                 ELSE
@@ -4420 +4423 @@
                       DO k = nzb+1, nzt+1
                          rrtm_h2ovmr(0,k) =  mol_mass_air_d_wv * q(k,j,i)
-                      ENDDO   
+                      ENDDO
                    ELSE
                       rrtm_h2ovmr(0,nzb+1:nzt+1) = 0.0_wp
@@ -4427 +4430 @@
 
 !
-!--             Avoid temperature/humidity jumps at the top of the LES domain by 
+!--             Avoid temperature/humidity jumps at the top of the LES domain by
 !--             linear interpolation from nzt+2 to nzt+7
                 DO k = nzt+2, nzt+7
@@ -4463 +4466 @@
                       rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *              &
                                           (rrtm_plev(0,k) - rrtm_plev(0,k+1))  &
-                                          * 100.0_wp / g 
+                                          * 100.0_wp / g
 
                       IF ( rrtm_cliqwp(0,k) > 0.0_wp )  THEN
@@ -4473 +4476 @@
                          IF ( bulk_cloud_model )  THEN
 !
-!--                         Calculete effective droplet radius. In case of using 
-!--                         cloud_scheme = 'morrison' and a non reasonable number 
-!--                         of cloud droplets the inital aerosol number  
+!--                         Calculete effective droplet radius. In case of using
+!--                         cloud_scheme = 'morrison' and a non reasonable number
+!--                         of cloud droplets the inital aerosol number
 !--                         concentration is considered.
                             IF ( microphysics_morrison )  THEN
@@ -4482 +4485 @@
                                ELSE
                                   nc_rad = na_init
-                               ENDIF 
+                               ENDIF
                             ELSE
                                nc_rad = nc_const
-                            ENDIF  
+                            ENDIF
 
                             rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)     &
@@ -4525 +4528 @@
 
 !
-!--             Write surface emissivity and surface temperature at current 
+!--             Write surface emissivity and surface temperature at current
 !--             surface element on RRTMG-shaped array.
 !--             Please note, as RRTMG is a single column model, surface attributes
-!--             are only obtained from horizontally aligned surfaces (for 
-!--             simplicity). Taking surface attributes from horizontal and 
-!--             vertical walls would lead to multiple solutions.  
+!--             are only obtained from horizontally aligned surfaces (for
+!--             simplicity). Taking surface attributes from horizontal and
+!--             vertical walls would lead to multiple solutions.
 !--             Moreover, for natural- and urban-type surfaces, several surface
-!--             classes can exist at a surface element next to each other. 
-!--             To obtain bulk parameters, apply a weighted average for these 
-!--             surfaces. 
+!--             classes can exist at a surface element next to each other.
+!--             To obtain bulk parameters, apply a weighted average for these
+!--             surfaces.
                 DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
                    rrtm_emis = surf_lsm_h%frac(m,ind_veg_wall)  *              &
                                surf_lsm_h%emissivity(m,ind_veg_wall)  +        &
                                surf_lsm_h%frac(m,ind_pav_green) *              &
-                               surf_lsm_h%emissivity(m,ind_pav_green) +        & 
+                               surf_lsm_h%emissivity(m,ind_pav_green) +        &
                                surf_lsm_h%frac(m,ind_wat_win)   *              &
                                surf_lsm_h%emissivity(m,ind_wat_win)
                    rrtm_tsfc = surf_lsm_h%pt_surface(m) * exner(nzb)
-                ENDDO             
+                ENDDO
                 DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
                    rrtm_emis = surf_usm_h%frac(m,ind_veg_wall)  *              &
                                surf_usm_h%emissivity(m,ind_veg_wall)  +        &
                                surf_usm_h%frac(m,ind_pav_green) *              &
-                               surf_usm_h%emissivity(m,ind_pav_green) +        & 
+                               surf_usm_h%emissivity(m,ind_pav_green) +        &
                                surf_usm_h%frac(m,ind_wat_win)   *              &
                                surf_usm_h%emissivity(m,ind_wat_win)
@@ -4562 +4565 @@
 !--                which are allocated on the exact size as the rrtmg_lw is called.
 !--                As one dimesion is allocated with zero size, compiler complains
-!--                that rank of the array does not match that of the 
-!--                assumed-shaped arguments in the RRTMG library. In order to 
-!--                avoid this, write to dummy arguments and give pass the entire 
-!--                dummy array. Seems to be the only existing work-around.  
+!--                that rank of the array does not match that of the
+!--                assumed-shaped arguments in the RRTMG library. In order to
+!--                avoid this, write to dummy arguments and give pass the entire
+!--                dummy array. Seems to be the only existing work-around.
                    ALLOCATE( rrtm_lw_taucld_dum(1:nbndlw+1,0:0,k_topo+1:nzt_rad+1) )
                    ALLOCATE( rrtm_lw_tauaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1) )
@@ -4574 +4577 @@
                                rrtm_lw_tauaer(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1)
 
-                   CALL rrtmg_lw( 1,                                           &                                        
+                   CALL rrtmg_lw( 1,                                           &
                                   nzt_rad-k_topo,                              &
                                   rrtm_icld,                                   &
@@ -4601 +4604 @@
                                   rrtm_cicewp(:,k_topo+1:nzt_rad+1),           &
                                   rrtm_cliqwp(:,k_topo+1:nzt_rad+1),           &
-                                  rrtm_reice(:,k_topo+1:nzt_rad+1),            & 
+                                  rrtm_reice(:,k_topo+1:nzt_rad+1),            &
                                   rrtm_reliq(:,k_topo+1:nzt_rad+1),            &
                                   rrtm_lw_tauaer_dum,                          &
@@ -4630 +4633 @@
 
 !
-!--                Save surface radiative fluxes and change in LW heating rate 
+!--                Save surface radiative fluxes and change in LW heating rate
 !--                onto respective surface elements
 !--                Horizontal surfaces
@@ -4638 +4641 @@
                       surf_lsm_h%rad_lw_out(m)          = rrtm_lwuflx(0,k_topo)
                       surf_lsm_h%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k_topo)
-                   ENDDO             
+                   ENDDO
                    DO  m = surf_usm_h%start_index(j,i),                        &
                            surf_usm_h%end_index(j,i)
@@ -4644 +4647 @@
                       surf_usm_h%rad_lw_out(m)          = rrtm_lwuflx(0,k_topo)
                       surf_usm_h%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k_topo)
-                   ENDDO 
-!
-!--                Vertical surfaces. Fluxes are obtain at vertical level of the 
+                   ENDDO
+!
+!--                Vertical surfaces. Fluxes are obtain at vertical level of the
 !--                respective surface element
                    DO  l = 0, 3
@@ -4655 +4658 @@
                          surf_lsm_v(l)%rad_lw_out(m)          = rrtm_lwuflx(0,k)
                          surf_lsm_v(l)%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k)
-                      ENDDO             
+                      ENDDO
                       DO  m = surf_usm_v(l)%start_index(j,i),                  &
                               surf_usm_v(l)%end_index(j,i)
@@ -4662 +4665 @@
                          surf_usm_v(l)%rad_lw_out(m)          = rrtm_lwuflx(0,k)
                          surf_usm_v(l)%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k)
-                      ENDDO 
+                      ENDDO
                    ENDDO
 
@@ -4669 +4672 @@
                 IF ( sw_radiation .AND. sun_up )  THEN
 !
-!--                Get albedo for direct/diffusive long/shortwave radiation at 
-!--                current (y,x)-location from surface variables. 
-!--                Only obtain it from horizontal surfaces, as RRTMG is a single 
+!--                Get albedo for direct/diffusive long/shortwave radiation at
+!--                current (y,x)-location from surface variables.
+!--                Only obtain it from horizontal surfaces, as RRTMG is a single
 !--                column model
-!--                (Please note, only one loop will entered, controlled by 
+!--                (Please note, only one loop will entered, controlled by
 !--                start-end index.)
                    DO  m = surf_lsm_h%start_index(j,i),                        &
                            surf_lsm_h%end_index(j,i)
                       rrtm_asdir(1)  = SUM( surf_lsm_h%frac(m,:) *             &
-                                            surf_lsm_h%rrtm_asdir(:,m) )
+                                            surf_lsm_h%rrtm_asdir(m,:) )
                       rrtm_asdif(1)  = SUM( surf_lsm_h%frac(m,:) *             &
-                                            surf_lsm_h%rrtm_asdif(:,m) )
+                                            surf_lsm_h%rrtm_asdif(m,:) )
                       rrtm_aldir(1)  = SUM( surf_lsm_h%frac(m,:) *             &
-                                            surf_lsm_h%rrtm_aldir(:,m) )
+                                            surf_lsm_h%rrtm_aldir(m,:) )
                       rrtm_aldif(1)  = SUM( surf_lsm_h%frac(m,:) *             &
-                                            surf_lsm_h%rrtm_aldif(:,m) )
-                   ENDDO             
+                                            surf_lsm_h%rrtm_aldif(m,:) )
+                   ENDDO
                    DO  m = surf_usm_h%start_index(j,i),                        &
                            surf_usm_h%end_index(j,i)
                       rrtm_asdir(1)  = SUM( surf_usm_h%frac(m,:) *             &
-                                            surf_usm_h%rrtm_asdir(:,m) )
+                                            surf_usm_h%rrtm_asdir(m,:) )
                       rrtm_asdif(1)  = SUM( surf_usm_h%frac(m,:) *             &
-                                            surf_usm_h%rrtm_asdif(:,m) )
+                                            surf_usm_h%rrtm_asdif(m,:) )
                       rrtm_aldir(1)  = SUM( surf_usm_h%frac(m,:) *             &
-                                            surf_usm_h%rrtm_aldir(:,m) )
+                                            surf_usm_h%rrtm_aldir(m,:) )
                       rrtm_aldif(1)  = SUM( surf_usm_h%frac(m,:) *             &
-                                            surf_usm_h%rrtm_aldif(:,m) )
+                                            surf_usm_h%rrtm_aldif(m,:) )
                    ENDDO
 !
-!--                Due to technical reasons, copy optical depths and other 
-!--                to dummy arguments which are allocated on the exact size as the 
+!--                Due to technical reasons, copy optical depths and other
+!--                to dummy arguments which are allocated on the exact size as the
 !--                rrtmg_sw is called.
 !--                As one dimesion is allocated with zero size, compiler complains
-!--                that rank of the array does not match that of the 
-!--                assumed-shaped arguments in the RRTMG library. In order to 
-!--                avoid this, write to dummy arguments and give pass the entire 
-!--                dummy array. Seems to be the only existing work-around.  
+!--                that rank of the array does not match that of the
+!--                assumed-shaped arguments in the RRTMG library. In order to
+!--                avoid this, write to dummy arguments and give pass the entire
+!--                dummy array. Seems to be the only existing work-around.
                    ALLOCATE( rrtm_sw_taucld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
                    ALLOCATE( rrtm_sw_ssacld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
@@ -4714 +4717 @@
                    ALLOCATE( rrtm_sw_asmaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
                    ALLOCATE( rrtm_sw_ecaer_dum(0:0,k_topo+1:nzt_rad+1,1:naerec+1)  )
-     
+
                    rrtm_sw_taucld_dum = rrtm_sw_taucld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
                    rrtm_sw_ssacld_dum = rrtm_sw_ssacld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
@@ -4733 +4736 @@
                                   rrtm_tlev(:,k_topo+1:nzt_rad+2),             &
                                   rrtm_tsfc,                                   &
-                                  rrtm_h2ovmr(:,k_topo+1:nzt_rad+1),           &                               
-                                  rrtm_o3vmr(:,k_topo+1:nzt_rad+1),            &        
+                                  rrtm_h2ovmr(:,k_topo+1:nzt_rad+1),           &
+                                  rrtm_o3vmr(:,k_topo+1:nzt_rad+1),            &
                                   rrtm_co2vmr(:,k_topo+1:nzt_rad+1),           &
                                   rrtm_ch4vmr(:,k_topo+1:nzt_rad+1),           &
                                   rrtm_n2ovmr(:,k_topo+1:nzt_rad+1),           &
                                   rrtm_o2vmr(:,k_topo+1:nzt_rad+1),            &
-                                  rrtm_asdir,                                  & 
+                                  rrtm_asdir,                                  &
                                   rrtm_asdif,                                  &
                                   rrtm_aldir,                                  &
@@ -4763 +4766 @@
                                   rrtm_sw_asmaer_dum,                          &
                                   rrtm_sw_ecaer_dum,                           &
-                                  rrtm_swuflx(:,k_topo:nzt_rad+1),             &  
-                                  rrtm_swdflx(:,k_topo:nzt_rad+1),             & 
-                                  rrtm_swhr(:,k_topo+1:nzt_rad+1),             & 
-                                  rrtm_swuflxc(:,k_topo:nzt_rad+1),            & 
+                                  rrtm_swuflx(:,k_topo:nzt_rad+1),             &
+                                  rrtm_swdflx(:,k_topo:nzt_rad+1),             &
+                                  rrtm_swhr(:,k_topo+1:nzt_rad+1),             &
+                                  rrtm_swuflxc(:,k_topo:nzt_rad+1),            &
                                   rrtm_swdflxc(:,k_topo:nzt_rad+1),            &
                                   rrtm_swhrc(:,k_topo+1:nzt_rad+1),            &
@@ -4800 +4803 @@
                       surf_lsm_h%rad_sw_in(m)     = rrtm_swdflx(0,k_topo)
                       surf_lsm_h%rad_sw_out(m)    = rrtm_swuflx(0,k_topo)
-                   ENDDO             
+                   ENDDO
                    DO  m = surf_usm_h%start_index(j,i),                        &
                            surf_usm_h%end_index(j,i)
                       surf_usm_h%rad_sw_in(m)     = rrtm_swdflx(0,k_topo)
                       surf_usm_h%rad_sw_out(m)    = rrtm_swuflx(0,k_topo)
-                   ENDDO 
-!
-!--                Vertical surfaces. Fluxes are obtain at respective vertical 
+                   ENDDO
+!
+!--                Vertical surfaces. Fluxes are obtain at respective vertical
 !--                level of the surface element
                    DO  l = 0, 3
@@ -4815 +4818 @@
                          surf_lsm_v(l)%rad_sw_in(m)  = rrtm_swdflx(0,k)
                          surf_lsm_v(l)%rad_sw_out(m) = rrtm_swuflx(0,k)
-                      ENDDO             
+                      ENDDO
                       DO  m = surf_usm_v(l)%start_index(j,i),                  &
                               surf_usm_v(l)%end_index(j,i)
@@ -4821 +4824 @@
                          surf_usm_v(l)%rad_sw_in(m)  = rrtm_swdflx(0,k)
                          surf_usm_v(l)%rad_sw_out(m) = rrtm_swuflx(0,k)
-                      ENDDO 
+                      ENDDO
                    ENDDO
 !
@@ -4829 +4832 @@
                    rad_sw_out = 0.0_wp
 !--             !!!!!!!! ATTENSION !!!!!!!!!!!!!!!
-!--             Surface radiative fluxes should be also set to zero here                  
+!--             Surface radiative fluxes should be also set to zero here
 !--                Save surface radiative fluxes onto respective surface elements
 !--                Horizontal surfaces
@@ -4836 +4839 @@
                       surf_lsm_h%rad_sw_in(m)     = 0.0_wp
                       surf_lsm_h%rad_sw_out(m)    = 0.0_wp
-                   ENDDO             
+                   ENDDO
                    DO  m = surf_usm_h%start_index(j,i),                        &
                            surf_usm_h%end_index(j,i)
                       surf_usm_h%rad_sw_in(m)     = 0.0_wp
                       surf_usm_h%rad_sw_out(m)    = 0.0_wp
-                   ENDDO 
-!
-!--                Vertical surfaces. Fluxes are obtain at respective vertical 
+                   ENDDO
+!
+!--                Vertical surfaces. Fluxes are obtain at respective vertical
 !--                level of the surface element
                    DO  l = 0, 3
@@ -4851 +4854 @@
                          surf_lsm_v(l)%rad_sw_in(m)  = 0.0_wp
                          surf_lsm_v(l)%rad_sw_out(m) = 0.0_wp
-                      ENDDO             
+                      ENDDO
                       DO  m = surf_usm_v(l)%start_index(j,i),                  &
                               surf_usm_v(l)%end_index(j,i)
@@ -4857 +4860 @@
                          surf_usm_v(l)%rad_sw_in(m)  = 0.0_wp
                          surf_usm_v(l)%rad_sw_out(m) = 0.0_wp
-                      ENDDO 
+                      ENDDO
                    ENDDO
                 ENDIF
@@ -4868 +4871 @@
 !--    Finally, calculate surface net radiation for surface elements.
        IF (  .NOT.  radiation_interactions  ) THEN
-!--       First, for horizontal surfaces    
+!--       First, for horizontal surfaces
           DO  m = 1, surf_lsm_h%ns
              surf_lsm_h%rad_net(m) = surf_lsm_h%rad_sw_in(m)                   &
@@ -4882 +4885 @@
           ENDDO
 !
-!--       Vertical surfaces. 
+!--       Vertical surfaces.
 !--       Todo: weight with azimuth and zenith angle according to their orientation!
-          DO  l = 0, 3     
+          DO  l = 0, 3
              DO  m = 1, surf_lsm_v(l)%ns
                 surf_lsm_v(l)%rad_net(m) = surf_lsm_v(l)%rad_sw_in(m)          &
@@ -4907 +4910 @@
 
        CALL exchange_horiz( rad_sw_in,  nbgp )
-       CALL exchange_horiz( rad_sw_out, nbgp ) 
+       CALL exchange_horiz( rad_sw_out, nbgp )
        CALL exchange_horiz( rad_sw_hr,    nbgp )
        CALL exchange_horiz( rad_sw_cs_hr, nbgp )
@@ -4921 +4924 @@
 !> Calculate the cosine of the zenith angle (variable is called zenith)
 !------------------------------------------------------------------------------!
-    SUBROUTINE calc_zenith( day_of_year, second_of_day ) 
+    SUBROUTINE calc_zenith( day_of_year, second_of_day )
 
        USE palm_date_time_mod,                                                 &
@@ -4973 +4976 @@
 ! Description:
 ! ------------
-!> Calculates surface albedo components based on Briegleb (1992) and 
+!> Calculates surface albedo components based on Briegleb (1992) and
 !> Briegleb et al. (1986)
 !------------------------------------------------------------------------------!
@@ -4990 +4993 @@
 !
 !--           Loop over surface elements
-              DO  ind_type = 0, SIZE( surf%albedo_type, 1 ) - 1
-           
+              DO  ind_type = 0, SIZE( surf%albedo_type, 2 ) - 1
+
 !
 !--              Ocean
@@ -5081 +5084 @@
            surf%rrtm_asdif = albedo_urb
            surf%rrtm_aldir = albedo_urb
-           surf%rrtm_aldif = albedo_urb  
+           surf%rrtm_aldif = albedo_urb
 !
 !--     Darkness
@@ -5150 +5153 @@
           DEALLOCATE ( rrtm_sw_tauaer )
           DEALLOCATE ( rrtm_sw_ssaaer )
-          DEALLOCATE ( rrtm_sw_asmaer ) 
-          DEALLOCATE ( rrtm_sw_ecaer )   
- 
+          DEALLOCATE ( rrtm_sw_asmaer )
+          DEALLOCATE ( rrtm_sw_ecaer )
+
           DEALLOCATE ( rrtm_swdflx  )
           DEALLOCATE ( rrtm_swdflxc )
@@ -5236 +5239 @@
 
 !
-!--    Calculate pressure levels on zu and zw grid. Sounding data is added at 
-!--    top of the LES domain. This routine does not consider horizontal or 
+!--    Calculate pressure levels on zu and zw grid. Sounding data is added at
+!--    top of the LES domain. This routine does not consider horizontal or
 !--    vertical variability of pressure and temperature
        ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1)   )
@@ -5263 +5266 @@
 
 !
-!--    Calculate temperature/humidity levels at top of the LES domain. 
-!--    Currently, the temperature is taken from sounding data (might lead to a 
-!--    temperature jump at interface. To do: Humidity is currently not 
+!--    Calculate temperature/humidity levels at top of the LES domain.
+!--    Currently, the temperature is taken from sounding data (might lead to a
+!--    temperature jump at interface. To do: Humidity is currently not
 !--    calculated above the LES domain.
        ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1)   )
@@ -5299 +5302 @@
        ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
        ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
-       ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) ) 
-       ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )    
+       ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
+       ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )
 
 !
@@ -5332 +5335 @@
        rrtm_swdflx  = 0.0_wp
        rrtm_swuflx  = 0.0_wp
-       rrtm_swhr    = 0.0_wp  
+       rrtm_swhr    = 0.0_wp
        rrtm_swuflxc = 0.0_wp
        rrtm_swdflxc = 0.0_wp
@@ -5348 +5351 @@
        rrtm_lwdflx  = 0.0_wp
        rrtm_lwuflx  = 0.0_wp
-       rrtm_lwhr    = 0.0_wp  
+       rrtm_lwhr    = 0.0_wp
        rrtm_lwuflxc = 0.0_wp
        rrtm_lwdflxc = 0.0_wp
@@ -5367 +5370 @@
 !> Read trace gas data from file and convert into trace gas paths / volume
 !> mixing ratios. If a user-defined input file is provided it needs to follow
-!> the convections used in RRTMG (see respective netCDF files shipped with 
+!> the convections used in RRTMG (see respective netCDF files shipped with
 !> RRTMG)
 !------------------------------------------------------------------------------!
@@ -5422 +5425 @@
           DEALLOCATE ( rrtm_cfc22vmr )
           DEALLOCATE ( rrtm_ccl4vmr  )
-          DEALLOCATE ( rrtm_h2ovmr  )     
+          DEALLOCATE ( rrtm_h2ovmr  )
        ENDIF
 
@@ -5446 +5449 @@
        nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
        CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
-       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np) 
+       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np)
        CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
 
        nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
        CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
-       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs ) 
+       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs )
        CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
-   
-
-!
-!--    Allocate pressure, and trace gas arrays     
+
+
+!
+!--    Allocate pressure, and trace gas arrays
        ALLOCATE( p_mls(1:np) )
-       ALLOCATE( trace_mls(1:num_trace_gases,1:np) ) 
-       ALLOCATE( trace_mls_tmp(1:nabs,1:np) ) 
+       ALLOCATE( trace_mls(1:num_trace_gases,1:np) )
+       ALLOCATE( trace_mls_tmp(1:nabs,1:np) )
 
 
@@ -5482 +5485 @@
 !
 !--       Replace missing values by zero
-          WHERE ( trace_mls(n,:) > 2.0_wp )  
+          WHERE ( trace_mls(n,:) > 2.0_wp )
              trace_mls(n,:) = 0.0_wp
           END WHERE
@@ -5497 +5500 @@
        ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
 
-       rrtm_play_tmp(1:nzt_rad)   = rrtm_play(0,1:nzt_rad) 
+       rrtm_play_tmp(1:nzt_rad)   = rrtm_play(0,1:nzt_rad)
        rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
        rrtm_play_tmp(nzt_rad+1)   = rrtm_plev(0,nzt_rad+1) * 0.5_wp
        rrtm_plev_tmp(nzt_rad+2)   = MIN( 1.0E-4_wp, 0.25_wp                    &
                                          * rrtm_plev(0,nzt_rad+1) )
- 
+
 !
 !--    Calculate trace gas path (zero at surface) with interpolation to the
@@ -5509 +5512 @@
 
        trace_mls_path(nzb+1,:) = 0.0_wp
-       
+
        DO k = nzb+2, nzt_rad+2
           DO m = 1, num_trace_gases
@@ -5516 +5519 @@
 !
 !--          When the pressure level is higher than the trace gas pressure
-!--          level, assume that 
-             IF ( rrtm_plev_tmp(k-1) > p_mls(1) )  THEN             
-                
+!--          level, assume that
+             IF ( rrtm_plev_tmp(k-1) > p_mls(1) )  THEN
+
                 trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1)     &
                                       * ( rrtm_plev_tmp(k-1)                   &
@@ -5526 +5529 @@
 
 !
-!--          Integrate for each sounding level from the contributing p_mls 
+!--          Integrate for each sounding level from the contributing p_mls
 !--          levels
              DO n = 2, np
@@ -5553 +5556 @@
              ENDDO
 
-             IF ( rrtm_plev_tmp(k) < p_mls(np) )  THEN 
+             IF ( rrtm_plev_tmp(k) < p_mls(np) )  THEN
                 trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np)    &
                                       * ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
                                           - rrtm_plev_tmp(k)                   &
-                                        ) / g 
-             ENDIF  
+                                        ) / g
+             ENDIF
           ENDDO
        ENDDO
@@ -5617 +5620 @@
 
                 rrtm_h2ovmr(0,:) = trace_path_tmp(:)
-                
+
              CASE DEFAULT
 
@@ -5772 +5775 @@
      INTEGER(iwp)                      ::  pc_box_dimshift    !< transform for best accuracy
      INTEGER(iwp), DIMENSION(0:3)      ::  reorder = (/ 1, 0, 3, 2 /)
-                                           
+
      REAL(wp), DIMENSION(3,3)          ::  mrot               !< grid rotation matrix (zyx)
      REAL(wp), DIMENSION(3,0:nsurf_type)::  vnorm             !< face direction normal vectors (zyx)
@@ -5797 +5800 @@
      REAL(wp)                          ::  area_surf          !< total area of surfaces in all processor
      REAL(wp)                          ::  area_hor           !< total horizontal area of domain in all processor
-#if defined( __parallel )     
+#if defined( __parallel )
      REAL(wp), DIMENSION(1:7)          ::  combine_allreduce   !< dummy array used to combine several MPI_ALLREDUCE calls
      REAL(wp), DIMENSION(1:7)          ::  combine_allreduce_l !< dummy array used to combine several MPI_ALLREDUCE calls
@@ -5850 +5853 @@
      ENDIF
 !
-!--  Split downwelling shortwave radiation into a diffuse and a direct part. 
+!--  Split downwelling shortwave radiation into a diffuse and a direct part.
 !--  Note, if radiation scheme is RRTMG or diffuse radiation is externally
 !--  prescribed, this is not required. Please note, in case of external
-!--  radiation, the clear-sky model is applied during spinup, so that 
-!--  radiation need to be split also in this case. 
+!--  radiation, the clear-sky model is applied during spinup, so that
+!--  radiation need to be split also in this case.
      IF ( radiation_scheme == 'constant'   .OR.                                &
           radiation_scheme == 'clear-sky'  .OR.                                &
@@ -6141 +6144 @@
      IF ( .NOT.  surface_reflections )  THEN
 !
-!--     Set nrefsteps to 0 to disable reflections        
+!--     Set nrefsteps to 0 to disable reflections
         nrefsteps = 0
         surfoutsl = albedo_surf * surfins
@@ -6549 +6552 @@
      combine_allreduce_l(6) = emiss_sum_surfl
      combine_allreduce_l(7) = area_surfl
-     
+
      CALL MPI_ALLREDUCE( combine_allreduce_l,                                  &
                          combine_allreduce,                                    &
@@ -6557 +6560 @@
                          comm2d,                                               &
                          ierr )
-     
+
      pinsw          = combine_allreduce(1)
      pinlw          = combine_allreduce(2)
@@ -6580 +6583 @@
      IF ( area_surf /= 0.0_wp )  emissivity_urb = emiss_sum_surf / area_surf
 !
-!--  Temporally comment out calculation of effective radiative temperature. 
+!--  Temporally comment out calculation of effective radiative temperature.
 !--  See below for more explanation.
 !--  (3) temperature
 !--   first we calculate an effective horizontal area to account for
 !--   the effect of vertical surfaces (which contributes to LW emission)
-!--   We simply use the ratio of the total LW to the incoming LW flux 
+!--   We simply use the ratio of the total LW to the incoming LW flux
       area_hor = pinlw / rad_lw_in_diff(nyn,nxl)
       t_rad_urb = ( ( pemitlw - pabslw + emissivity_urb * pinlw ) / &
-           (emissivity_urb * sigma_sb * area_hor) )**0.25_wp     
+           (emissivity_urb * sigma_sb * area_hor) )**0.25_wp
 
      IF ( debug_output_timestep )  CALL debug_message( 'radiation_interaction', 'end' )
@@ -6675 +6678 @@
 !> It follows Boland, Ridley & Brown (2008)
 !------------------------------------------------------------------------------!
-    SUBROUTINE calc_diffusion_radiation 
+    SUBROUTINE calc_diffusion_radiation
 
        USE palm_date_time_mod,                                                 &
@@ -6706 +6709 @@
                           0.000719_wp * cos(2.0_wp * year_angle) +             &
                           0.000077_wp * sin(2.0_wp * year_angle))
-        
-!--    
+
+!--
 !--     Under a very low angle, we keep extraterestrial radiation at
 !--     the last small value, therefore the clearness index will be pushed
-!--     towards 0 while keeping full continuity.    
+!--     towards 0 while keeping full continuity.
         IF ( cos_zenith <= lowest_solarUp )  THEN
             corrected_solarUp = lowest_solarUp
@@ -6716 +6719 @@
             corrected_solarUp = cos_zenith
         ENDIF
-        
+
         horizontalETR = etr * corrected_solarUp
-        
+
         DO i = nxl, nxr
             DO j = nys, nyn
@@ -6728 +6731 @@
             ENDDO
         ENDDO
-        
+
     END SUBROUTINE calc_diffusion_radiation
 
@@ -6826 +6829 @@
 
  END SUBROUTINE radiation_interaction
-    
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -6934 +6937 @@
 #endif
 !
-!--    Stretching (non-uniform grid spacing) is not considered in the radiation 
-!--    model. Therefore, vertical stretching has to be applied above the area 
+!--    Stretching (non-uniform grid spacing) is not considered in the radiation
+!--    model. Therefore, vertical stretching has to be applied above the area
 !--    where the parts of the radiation model which assume constant grid spacing
-!--    are active. ABS (...) is required because the default value of 
-!--    dz_stretch_level_start is -9999999.9_wp (negative). 
+!--    are active. ABS (...) is required because the default value of
+!--    dz_stretch_level_start is -9999999.9_wp (negative).
        IF ( ABS( dz_stretch_level_start(1) ) <= zw(nz_urban_t) ) THEN
           WRITE( message_string, * ) 'The lowest level where vertical ',       &
@@ -6944 +6947 @@
                                      'greater than ', zw(nz_urban_t)
           CALL message( 'radiation_interaction_init', 'PA0496', 1, 2, 0, 6, 0 )
-       ENDIF 
+       ENDIF
 !
 !--    global number of urban and plant layers
@@ -7365 +7368 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE radiation_calc_svf
-    
+
         IMPLICIT NONE
-        
+
         INTEGER(iwp)                                  :: i, j, k, d, ip, jp
         INTEGER(iwp)                                  :: isvf, ksvf, icsf, kcsf, npcsfl, isvf_surflt, imrt, imrtf, ipcgb
@@ -7411 +7414 @@
         INTEGER(kind=MPI_ADDRESS_KIND)                :: size_lad_rma
 #endif
-!    
+!
         INTEGER(iwp), DIMENSION(0:svfnorm_report_num) :: svfnorm_counts
 
@@ -7440 +7443 @@
         ray_skip_maxdist = 0
         ray_skip_minval = 0
-        
+
 !--     initialize temporary terrain and plant canopy height arrays (global 2D array!)
         ALLOCATE( nzterr(0:(nx+1)*(ny+1)-1) )
@@ -7482 +7485 @@
             IF ( raytrace_mpi_rma )  THEN
                 ALLOCATE( lad_s_ray(maxboxesg) )
-                
+
                 ! set conditions for RMA communication
                 CALL MPI_Info_create(minfo, ierr)
@@ -7556 +7559 @@
                     FLUSH(9)
                 ENDIF
-                
+
             ELSE
                 ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nz_plant-1) )
@@ -7589 +7592 @@
         !--
         dsitrans(:,:) = 0._wp
-        
+
         DO isurflt = 1, nsurfl
 !--         determine face centers
@@ -7694 +7697 @@
                                  SUM(ztransp(itarg0:itarg0+lowest_free_ray-1) &
                                              * vffrac(itarg0:itarg0+lowest_free_ray-1))
- 
+
 !--            Save direct solar transparency
                j = MODULO(NINT(azmid/                                          &
@@ -7760 +7763 @@
                         CALL debug_message( debug_string, 'info' )
                      ENDIF
-                     
+
                      nsvfla = k
                   ENDIF
@@ -7785 +7788 @@
                      CYCLE
                   ENDIF
-                  
+
                   sd = surf(id, isurfs)
                   sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd),  &
@@ -7801 +7804 @@
                      CYCLE
                   ENDIF
-                 
+
                   difvf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction
                       * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) &  ! cosine of target normal and reverse direction
@@ -7845 +7848 @@
                         CALL debug_message( debug_string, 'info' )
                      ENDIF
-                     
+
                      nsvfla = k
                   ENDIF
@@ -8082 +8085 @@
 !--     deallocate temporary global arrays
         DEALLOCATE(nzterr)
-        
+
         IF ( plant_canopy )  THEN
 !--         finalize mpi_rma communication and deallocate temporary arrays
@@ -8202 +8205 @@
 !--         sort and merge csf for the last time, keeping the array size to minimum
             CALL merge_and_grow_csf(-1)
-            
+
 !--         aggregate csb among processors
 !--         allocate necessary arrays
@@ -8214 +8217 @@
             ALLOCATE( ipcsflt(0:numprocs-1) )
             ALLOCATE( dpcsflt(0:numprocs-1) )
-            
-!--         fill out arrays of csf values and 
+
+!--         fill out arrays of csf values and
 !--         arrays of number of elements and displacements
 !--         for particular precessors
@@ -8257 +8260 @@
                 dcsflt(jp) = d
             ENDDO
-            
+
 !--         deallocate temporary acsf array
 !--         DEALLOCATE(acsf) - ifort has a problem with deallocation of allocatable target
@@ -8267 +8270 @@
                 DEALLOCATE(acsf2)
             ENDIF
-                    
+
 #if defined( __parallel )
 !--         scatter and gather the number of elements to and from all processor
@@ -8309 +8312 @@
                 FLUSH(9)
             ENDIF
-            
+
 #else
             npcsfl = ncsfl
@@ -8370 +8373 @@
                 ENDDO
             ENDIF
-            
+
 !--         deallocation of temporary arrays
             IF ( npcbl > 0 )  DEALLOCATE( gridpcbl )
@@ -8379 +8382 @@
                CALL debug_message( debug_string, 'info' )
             ENDIF
-            
+
         ENDIF
 
@@ -8388 +8391 @@
 
         RETURN  !todo: remove
-        
+
 !        WRITE( message_string, * )  &
 !            'I/O error when processing shape view factors / ',  &
 !            'plant canopy sink factors / direct irradiance factors.'
 !        CALL message( 'init_urban_surface', 'PA0502', 2, 2, 0, 6, 0 )
-       
+
     END SUBROUTINE radiation_calc_svf
 
-    
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -8440 +8443 @@
         INTEGER(iwp), DIMENSION(3)             :: dimnext      !< next dimension increments along path
         INTEGER(iwp), DIMENSION(3)             :: dimdelta     !< dimension direction = +- 1
-        INTEGER(iwp)                           :: px, py       !< number of processors in x and y dir before 
+        INTEGER(iwp)                           :: px, py       !< number of processors in x and y dir before
                                                                !< the processor in the question
         INTEGER(iwp)                           :: ip           !< number of processor where gridbox reside
         INTEGER(iwp)                           :: ig           !< 1D index of gridbox in global 2D array
-        
+
         REAL(wp)                               :: eps = 1E-10_wp !< epsilon for value comparison
         REAL(wp)                               :: lad_s_target   !< recieved lad_s of particular grid box
@@ -8461 +8464 @@
             CALL merge_and_grow_csf(k)
         ENDIF
-        
+
         transparency = 1._wp
         ncsb = 0
@@ -8534 +8537 @@
             dimnextdist(seldim) = (dimnext(seldim) - .5_wp - src(seldim)) / uvect(seldim)
         ENDDO
-        
+
         IF ( plant_canopy )  THEN
 #if defined( __parallel )
@@ -8549 +8552 @@
                     ENDIF
                 ENDDO
-                
+
 !--             wait for all pending local requests complete
                 CALL MPI_Win_flush_local_all(win_lad, ierr)
@@ -8557 +8560 @@
                 ENDIF
                 CALL cpu_log( log_point_s(77), 'rad_rma_lad', 'stop' )
-                
+
             ENDIF
 #endif
@@ -8586 +8589 @@
 
                 transparency = transparency * (1._wp - cursink)
-                
+
             ENDDO
         ENDIF
-        
+
         visible = .TRUE.
 
     END SUBROUTINE raytrace
-    
- 
+
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -8640 +8643 @@
       INTEGER(iwp), DIMENSION(2)             ::  dimnext      !< next dimension increments along path
       INTEGER(iwp), DIMENSION(2)             ::  dimdelta     !< dimension direction = +- 1
-      INTEGER(iwp)                           ::  px, py       !< number of processors in x and y dir before 
+      INTEGER(iwp)                           ::  px, py       !< number of processors in x and y dir before
                                                               !< the processor in the question
       INTEGER(iwp)                           ::  ip           !< number of processor where gridbox reside
       INTEGER(iwp)                           ::  ig           !< 1D index of gridbox in global 2D array
       INTEGER(iwp)                           ::  maxboxes     !< max no of CSF created
-      INTEGER(iwp)                           ::  nly          !< maximum  plant canopy height 
+      INTEGER(iwp)                           ::  nly          !< maximum  plant canopy height
       INTEGER(iwp)                           ::  ntrack
-      
+
       INTEGER(iwp)                           ::  zb0
       INTEGER(iwp)                           ::  zb1
@@ -8661 +8664 @@
       INTEGER(MPI_ADDRESS_KIND)              ::  wdisp        !< RMA window displacement
 #endif
-      
+
       REAL(wp)                               ::  eps = 1E-10_wp !< epsilon for value comparison
       REAL(wp)                               ::  zbottom, ztop !< urban surface boundary in real numbers
@@ -8669 +8672 @@
       REAL(wp)                               ::  dxxyy         !< square of real horizontal distance
       REAL(wp)                               ::  curtrans      !< transparency of current PC box crossing
-      
-
-      
+
+
+
       yxorigin(:) = origin(2:3)
       transparency(:) = 1._wp !-- Pre-set the all rays to transparent before reducing
@@ -8816 +8819 @@
       IF ( plant_canopy )  THEN
 !--      Request LAD WHERE applicable
-!--      
+!--
 #if defined( __parallel )
          IF ( raytrace_mpi_rma )  THEN
@@ -8898 +8901 @@
             ENDIF
          ENDDO
-         
+
          DEALLOCATE( target_surfl )
-         
+
       ELSE
          itarget(:) = -1
@@ -8907 +8910 @@
       IF ( plant_canopy )  THEN
 !--      Skip the PCB around origin if requested (for MRT, the PCB might not be there)
-!--      
+!--
          IF ( skip_1st_pcb  .AND.  NINT(origin(1)) <= plantt_max )  THEN
             rt2_track_lad(NINT(origin(1), iwp), 1) = 0._wp
@@ -8913 +8916 @@
 
 !--      Assert that we have space allocated for CSFs
-!--      
+!--
          maxboxes = (ntrack + MAX(CEILING(origin(1)-.5_wp) - nz_urban_b,          &
                                   nz_urban_t - CEILING(origin(1)-.5_wp))) * nrays
@@ -8926 +8929 @@
 
 !--      Calculate transparencies and store new CSFs
-!--      
+!--
          zbottom = REAL(nz_urban_b, wp) - .5_wp
          ztop = REAL(plantt_max, wp) + .5_wp
 
 !--      Reverse direction of radiation (face->sky), only when calc_svf
-!--      
+!--
          IF ( calc_svf )  THEN
             DO  i = 1, ntrack ! for each column
@@ -9002 +9005 @@
 
 !--      Forward direction of radiation (sky->face), always
-!--      
+!--
          DO  i = ntrack, 1, -1 ! for each column backwards
             dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)-rt2_track_dist(i-1))**2
@@ -9107 +9110 @@
 
    END SUBROUTINE raytrace_2d
- 
+
 
 !------------------------------------------------------------------------------!
@@ -9132 +9135 @@
       ndsidir = 0
       sun_direction = .TRUE.
-      
+
 !
 !--   Process spinup time if configured
@@ -9176 +9179 @@
          CALL calc_zenith( day_of_year, second_of_day )
          IF ( cos_zenith > 0 )  THEN
-!--         
+!--
 !--         Identify solar direction vector (discretized number) 1)
             i = MODULO(NINT(ATAN2(sun_dir_lon, sun_dir_lat)               &
@@ -9210 +9213 @@
         IMPLICIT NONE
         INTEGER(iwp),   INTENT(in)  :: x, y, z, d, x2, y2, z2, d2
-      
+
         surface_facing = .FALSE.
 
@@ -9255 +9258 @@
 
         surface_facing = .TRUE.
-        
+
     END FUNCTION surface_facing
 
@@ -9270 +9273 @@
 
        IMPLICIT NONE
-       
+
        CHARACTER(rad_version_len)   :: rad_version_field
-       
+
        INTEGER(iwp)                 :: i
        INTEGER(iwp)                 :: ndsidir_from_file = 0
@@ -9296 +9299 @@
                                               'that has written the svf data ',&
                                               'and the one that will read it ',&
-                                              'is not allowed' 
+                                              'is not allowed'
                    CALL message( 'check_open', 'PA0491', 1, 2, 0, 6, 0 )
                 ENDIF
 
              ENDIF
-              
-!
-!--          Open binary file 
+
+!
+!--          Open binary file
              CALL check_open( 88 )
 
@@ -9315 +9318 @@
                  CALL message( 'radiation_read_svf', 'PA0482', 1, 2, 0, 6, 0 )
              ENDIF
-              
+
 !
 !--          read nsvfl, ncsfl, nsurfl, nmrtf
              READ ( 88 ) nsvfl, ncsfl, nsurfl_from_file, npcbl_from_file,      &
                          ndsidir_from_file, nmrtbl_from_file, nmrtf
-             
+
              IF ( nsvfl < 0  .OR.  ncsfl < 0 )  THEN
                  WRITE( message_string, * ) 'Wrong number of SVF or CSF'
@@ -9330 +9333 @@
                  IF ( debug_output )  CALL debug_message( debug_string, 'info' )
              ENDIF
-             
+
              IF ( nsurfl_from_file /= nsurfl )  THEN
                  WRITE( message_string, * ) 'nsurfl from SVF file does not ',  &
@@ -9337 +9340 @@
                  CALL message( 'radiation_read_svf', 'PA0490', 1, 2, 0, 6, 0 )
              ENDIF
-             
+
              IF ( npcbl_from_file /= npcbl )  THEN
                  WRITE( message_string, * ) 'npcbl from SVF file does not ',   &
@@ -9344 +9347 @@
                  CALL message( 'radiation_read_svf', 'PA0493', 1, 2, 0, 6, 0 )
              ENDIF
-             
+
              IF ( ndsidir_from_file /= ndsidir )  THEN
                  WRITE( message_string, * ) 'ndsidir from SVF file does not ', &
@@ -9360 +9363 @@
                  IF ( debug_output )  CALL debug_message( debug_string, 'info' )
              ENDIF
-              
-!
-!--          Arrays skyvf, skyvft, dsitrans and dsitransc are allready 
-!--          allocated in radiation_interaction_init and 
+
+!
+!--          Arrays skyvf, skyvft, dsitrans and dsitransc are allready
+!--          allocated in radiation_interaction_init and
 !--          radiation_presimulate_solar_pos
              IF ( nsurfl > 0 )  THEN
                 READ(88) skyvf
                 READ(88) skyvft
-                READ(88) dsitrans  
+                READ(88) dsitrans
              ENDIF
-              
+
              IF ( plant_canopy  .AND.  npcbl > 0 ) THEN
                 READ ( 88 )  dsitransc
              ENDIF
-             
-!
-!--          The allocation of svf, svfsurf, csf, csfsurf, mrtf, mrtft, and 
-!--          mrtfsurf happens in routine radiation_calc_svf which is not 
-!--          called if the program enters radiation_read_svf. Therefore 
+
+!
+!--          The allocation of svf, svfsurf, csf, csfsurf, mrtf, mrtft, and
+!--          mrtfsurf happens in routine radiation_calc_svf which is not
+!--          called if the program enters radiation_read_svf. Therefore
 !--          these arrays has to allocate in the following
              IF ( nsvfl > 0 )  THEN
@@ -9408 +9411 @@
                 READ(88) mrtfsurf
              ENDIF
-              
-!
-!--          Close binary file                 
+
+!
+!--          Close binary file
              CALL close_file( 88 )
-                
+
           ENDIF
 #if defined( __parallel )
@@ -9436 +9439 @@
 
        IMPLICIT NONE
-       
+
        INTEGER(iwp)        :: i
 
@@ -9445 +9448 @@
           IF ( i == io_group )  THEN
 !
-!--          Open binary file 
+!--          Open binary file
              CALL check_open( 89 )
 
@@ -9473 +9476 @@
              IF ( nmrtf > 0 )  THEN
                  WRITE ( 89 )  mrtf
-                 WRITE ( 89 )  mrtft                
+                 WRITE ( 89 )  mrtft
                  WRITE ( 89 )  mrtfsurf
              ENDIF
 !
-!--          Close binary file                 
+!--          Close binary file
              CALL close_file( 89 )
 
@@ -9627 +9630 @@
     END SUBROUTINE quicksort_csf
 
-    
+
 !------------------------------------------------------------------------------!
 !
@@ -9712 +9715 @@
     END SUBROUTINE merge_and_grow_csf
 
-    
+
 !-- quicksort.f -*-f90-*-
 !-- Author: t-nissie, adaptation J.Resler
@@ -9746 +9749 @@
         IF ( j+1 < last )  CALL quicksort_csf2(kpcsflt, pcsflt, j+1, last)
     END SUBROUTINE quicksort_csf2
-    
+
 
     PURE FUNCTION csf_lt2(item1, item2) result(res)
@@ -9781 +9784 @@
 !------------------------------------------------------------------------------!
 SUBROUTINE radiation_3d_data_averaging( mode, variable )
- 
+
 
     USE control_parameters
@@ -9791 +9794 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  mode    !< 
-    CHARACTER (LEN=*) :: variable !< 
+    CHARACTER (LEN=*) ::  mode    !<
+    CHARACTER (LEN=*) :: variable !<
 
     LOGICAL      ::  match_lsm !< flag indicating natural-type surface
     LOGICAL      ::  match_usm !< flag indicating urban-type surface
-    
-    INTEGER(iwp) ::  i !< 
-    INTEGER(iwp) ::  j !< 
-    INTEGER(iwp) ::  k !< 
+
+    INTEGER(iwp) ::  i !<
+    INTEGER(iwp) ::  j !<
+    INTEGER(iwp) ::  k !<
     INTEGER(iwp) ::  l, m !< index of current surface element
 
@@ -9835 +9838 @@
                 ENDIF
                 rad_net_av = 0.0_wp
-             
+
              CASE ( 'rad_lw_in*' )
                 IF ( .NOT. ALLOCATED( rad_lw_in_xy_av ) )  THEN
@@ -9841 +9844 @@
                 ENDIF
                 rad_lw_in_xy_av = 0.0_wp
-                
+
              CASE ( 'rad_lw_out*' )
                 IF ( .NOT. ALLOCATED( rad_lw_out_xy_av ) )  THEN
@@ -9847 +9850 @@
                 ENDIF
                 rad_lw_out_xy_av = 0.0_wp
-                
+
              CASE ( 'rad_sw_in*' )
                 IF ( .NOT. ALLOCATED( rad_sw_in_xy_av ) )  THEN
@@ -9853 +9856 @@
                 ENDIF
                 rad_sw_in_xy_av = 0.0_wp
-                
+
              CASE ( 'rad_sw_out*' )
                 IF ( .NOT. ALLOCATED( rad_sw_out_xy_av ) )  THEN
                    ALLOCATE( rad_sw_out_xy_av(nysg:nyng,nxlg:nxrg) )
                 ENDIF
-                rad_sw_out_xy_av = 0.0_wp                
+                rad_sw_out_xy_av = 0.0_wp
 
              CASE ( 'rad_lw_in' )
@@ -10097 +10100 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'rad_lw_out*' )
              IF ( ALLOCATED( rad_lw_out_xy_av ) ) THEN
@@ -10119 +10122 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'rad_sw_in*' )
              IF ( ALLOCATED( rad_sw_in_xy_av ) ) THEN
@@ -10141 +10144 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'rad_sw_out*' )
              IF ( ALLOCATED( rad_sw_out_xy_av ) ) THEN
@@ -10163 +10166 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'rad_lw_in' )
              IF ( ALLOCATED( rad_lw_in_av ) ) THEN
@@ -10418 +10421 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'rad_lw_in*' )
              IF ( ALLOCATED( rad_lw_in_xy_av ) ) THEN
@@ -10428 +10431 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'rad_lw_out*' )
              IF ( ALLOCATED( rad_lw_out_xy_av ) ) THEN
@@ -10438 +10441 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'rad_sw_in*' )
              IF ( ALLOCATED( rad_sw_in_xy_av ) ) THEN
@@ -10448 +10451 @@
                 ENDDO
              ENDIF
-             
+
           CASE ( 'rad_sw_out*' )
              IF ( ALLOCATED( rad_sw_out_xy_av ) ) THEN
@@ -10702 +10705 @@
 !------------------------------------------------------------------------------!
 SUBROUTINE radiation_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
-    
+
     IMPLICIT NONE
 
     CHARACTER (LEN=*), INTENT(IN)  ::  variable    !<
-    LOGICAL, INTENT(OUT)           ::  found       !< 
-    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !< 
-    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !< 
-    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !< 
+    LOGICAL, INTENT(OUT)           ::  found       !<
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
 
     CHARACTER (len=varnamelength)  :: var
@@ -10779 +10782 @@
  SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode,         &
                                       local_pf, two_d, nzb_do, nzt_do )
- 
+
     USE indices
 
@@ -10787 +10790 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  grid     !< 
-    CHARACTER (LEN=*) ::  mode     !< 
-    CHARACTER (LEN=*) ::  variable !< 
-
-    INTEGER(iwp) ::  av !< 
-    INTEGER(iwp) ::  i  !< 
-    INTEGER(iwp) ::  j  !< 
-    INTEGER(iwp) ::  k  !< 
+    CHARACTER (LEN=*) ::  grid     !<
+    CHARACTER (LEN=*) ::  mode     !<
+    CHARACTER (LEN=*) ::  variable !<
+
+    INTEGER(iwp) ::  av !<
+    INTEGER(iwp) ::  i  !<
+    INTEGER(iwp) ::  j  !<
+    INTEGER(iwp) ::  k  !<
     INTEGER(iwp) ::  m  !< index of surface element at grid point (j,i)
-    INTEGER(iwp) ::  nzb_do   !< 
-    INTEGER(iwp) ::  nzt_do   !< 
-
-    LOGICAL      ::  found !< 
+    INTEGER(iwp) ::  nzb_do   !<
+    INTEGER(iwp) ::  nzt_do   !<
+
+    LOGICAL      ::  found !<
     LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
 
     REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
 
-    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
+    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
 
     found = .TRUE.
@@ -10818 +10821 @@
 !--                Natural-type surfaces
                    DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)  
+                           surf_lsm_h%end_index(j,i)
                       local_pf(i,j,nzb+1) = surf_lsm_h%rad_net(m)
                    ENDDO
@@ -10824 +10827 @@
 !--                Urban-type surfaces
                    DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)  
+                           surf_usm_h%end_index(j,i)
                       local_pf(i,j,nzb+1) = surf_usm_h%rad_net(m)
                    ENDDO
@@ -10835 +10838 @@
              ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    local_pf(i,j,nzb+1) = rad_net_av(j,i)
                 ENDDO
@@ -10842 +10845 @@
           two_d = .TRUE.
           grid = 'zu1'
-          
+
        CASE ( 'rad_lw_in*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
@@ -10851 +10854 @@
 !--                Natural-type surfaces
                    DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)  
+                           surf_lsm_h%end_index(j,i)
                       local_pf(i,j,nzb+1) = surf_lsm_h%rad_lw_in(m)
                    ENDDO
@@ -10857 +10860 @@
 !--                Urban-type surfaces
                    DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)  
+                           surf_usm_h%end_index(j,i)
                       local_pf(i,j,nzb+1) = surf_usm_h%rad_lw_in(m)
                    ENDDO
@@ -10868 +10871 @@
              ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    local_pf(i,j,nzb+1) = rad_lw_in_xy_av(j,i)
                 ENDDO
@@ -10875 +10878 @@
           two_d = .TRUE.
           grid = 'zu1'
-          
+
        CASE ( 'rad_lw_out*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
@@ -10884 +10887 @@
 !--                Natural-type surfaces
                    DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)  
+                           surf_lsm_h%end_index(j,i)
                       local_pf(i,j,nzb+1) = surf_lsm_h%rad_lw_out(m)
                    ENDDO
@@ -10890 +10893 @@
 !--                Urban-type surfaces
                    DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)  
+                           surf_usm_h%end_index(j,i)
                       local_pf(i,j,nzb+1) = surf_usm_h%rad_lw_out(m)
                    ENDDO
@@ -10901 +10904 @@
              ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    local_pf(i,j,nzb+1) = rad_lw_out_xy_av(j,i)
                 ENDDO
@@ -10908 +10911 @@
           two_d = .TRUE.
           grid = 'zu1'
-          
+
        CASE ( 'rad_sw_in*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
@@ -10917 +10920 @@
 !--                Natural-type surfaces
                    DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)  
+                           surf_lsm_h%end_index(j,i)
                       local_pf(i,j,nzb+1) = surf_lsm_h%rad_sw_in(m)
                    ENDDO
@@ -10923 +10926 @@
 !--                Urban-type surfaces
                    DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)  
+                           surf_usm_h%end_index(j,i)
                       local_pf(i,j,nzb+1) = surf_usm_h%rad_sw_in(m)
                    ENDDO
@@ -10934 +10937 @@
              ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    local_pf(i,j,nzb+1) = rad_sw_in_xy_av(j,i)
                 ENDDO
@@ -10941 +10944 @@
           two_d = .TRUE.
           grid = 'zu1'
-          
+
        CASE ( 'rad_sw_out*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
@@ -10950 +10953 @@
 !--                Natural-type surfaces
                    DO  m = surf_lsm_h%start_index(j,i),                        &
-                           surf_lsm_h%end_index(j,i)  
+                           surf_lsm_h%end_index(j,i)
                       local_pf(i,j,nzb+1) = surf_lsm_h%rad_sw_out(m)
                    ENDDO
@@ -10956 +10959 @@
 !--                Urban-type surfaces
                    DO  m = surf_usm_h%start_index(j,i),                        &
-                           surf_usm_h%end_index(j,i)  
+                           surf_usm_h%end_index(j,i)
                       local_pf(i,j,nzb+1) = surf_usm_h%rad_sw_out(m)
                    ENDDO
@@ -10967 +10970 @@
              ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    local_pf(i,j,nzb+1) = rad_sw_out_xy_av(j,i)
                 ENDDO
@@ -10973 +10976 @@
           ENDIF
           two_d = .TRUE.
-          grid = 'zu1'          
-          
+          grid = 'zu1'
+
        CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' )
           IF ( av == 0 ) THEN
@@ -10990 +10993 @@
             ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_lw_in_av(k,j,i)
@@ -11014 +11017 @@
             ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_lw_out_av(k,j,i)
@@ -11020 +11023 @@
                 ENDDO
              ENDDO
-          ENDIF   
+          ENDIF
           IF ( mode == 'xy' )  grid = 'zu'
 
@@ -11038 +11041 @@
             ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
@@ -11062 +11065 @@
             ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
@@ -11086 +11089 @@
             ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_sw_in_av(k,j,i)
@@ -11110 +11113 @@
             ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    DO  k = nzb, nzt+1
                       local_pf(i,j,k) = rad_sw_out_av(k,j,i)
@@ -11134 +11137 @@
             ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
@@ -11158 +11161 @@
             ENDIF
              DO  i = nxl, nxr
-                DO  j = nys, nyn 
+                DO  j = nys, nyn
                    DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
@@ -11172 +11175 @@
 
     END SELECT
- 
+
  END SUBROUTINE radiation_data_output_2d
 
@@ -11183 +11186 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
- 
+
 
     USE indices
@@ -11192 +11195 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  variable !< 
+    CHARACTER (LEN=*) ::  variable !<
 
     INTEGER(iwp) ::  av          !<
@@ -11217 +11220 @@
        RETURN
     ENDIF
-    
+
     IF ( var(1:3) /= 'rad'  .AND.  var(1:3) /= 'rtm' )  THEN
        found = .FALSE.
@@ -11684 +11687 @@
             ENDIF
          ENDIF
-!         
+!
 !--   block of RTM output variables
 !--   variables are intended mainly for debugging and detailed analyse purposes
       CASE ( 'rtm_skyvf' )
-!      
+!
 !--      sky view factor
          DO isurf = dirstart(ids), dirend(ids)
@@ -11758 +11761 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE radiation_data_output_mask( av, variable, found, local_pf, mid )
- 
+
     USE control_parameters
-        
+
     USE indices
-    
+
     USE kinds
-    
+
 
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  variable   !< 
+    CHARACTER (LEN=*) ::  variable   !<
 
     CHARACTER(LEN=5) ::  grid        !< flag to distinquish between staggered grids
 
-    INTEGER(iwp) ::  av              !< 
-    INTEGER(iwp) ::  i               !< 
-    INTEGER(iwp) ::  j               !< 
-    INTEGER(iwp) ::  k               !< 
+    INTEGER(iwp) ::  av              !<
+    INTEGER(iwp) ::  i               !<
+    INTEGER(iwp) ::  j               !<
+    INTEGER(iwp) ::  k               !<
     INTEGER(iwp) ::  mid             !< masked output running index
     INTEGER(iwp) ::  topo_top_index  !< k index of highest horizontal surface
@@ -11785 +11788 @@
     REAL(wp),                                                                  &
        DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
-          local_pf   !< 
+          local_pf   !<
 
     REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which needs to be resorted for output
@@ -11918 +11921 @@
        WRITE ( 14 )  rad_net_av
     ENDIF
-    
+
     IF ( ALLOCATED( rad_lw_in_xy_av ) )  THEN
        CALL wrd_write_string( 'rad_lw_in_xy_av' )
        WRITE ( 14 )  rad_lw_in_xy_av
     ENDIF
-    
+
     IF ( ALLOCATED( rad_lw_out_xy_av ) )  THEN
        CALL wrd_write_string( 'rad_lw_out_xy_av' )
        WRITE ( 14 )  rad_lw_out_xy_av
     ENDIF
-    
+
     IF ( ALLOCATED( rad_sw_in_xy_av ) )  THEN
        CALL wrd_write_string( 'rad_sw_in_xy_av' )
        WRITE ( 14 )  rad_sw_in_xy_av
     ENDIF
-    
+
     IF ( ALLOCATED( rad_sw_out_xy_av ) )  THEN
        CALL wrd_write_string( 'rad_sw_out_xy_av' )
@@ -12030 +12033 @@
                                 nxr_on_file, nynf, nync, nyn_on_file, nysf,    &
                                 nysc, nys_on_file, tmp_2d, tmp_3d, found )
- 
+
 
     USE control_parameters
-        
+
     USE indices
-    
+
     USE kinds
-    
+
     USE pegrid
 
@@ -12043 +12046 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  k               !< 
-    INTEGER(iwp) ::  nxlc            !< 
-    INTEGER(iwp) ::  nxlf            !< 
-    INTEGER(iwp) ::  nxl_on_file     !< 
-    INTEGER(iwp) ::  nxrc            !< 
-    INTEGER(iwp) ::  nxrf            !< 
-    INTEGER(iwp) ::  nxr_on_file     !< 
-    INTEGER(iwp) ::  nync            !< 
-    INTEGER(iwp) ::  nynf            !< 
-    INTEGER(iwp) ::  nyn_on_file     !< 
-    INTEGER(iwp) ::  nysc            !< 
-    INTEGER(iwp) ::  nysf            !< 
-    INTEGER(iwp) ::  nys_on_file     !< 
+    INTEGER(iwp) ::  k               !<
+    INTEGER(iwp) ::  nxlc            !<
+    INTEGER(iwp) ::  nxlf            !<
+    INTEGER(iwp) ::  nxl_on_file     !<
+    INTEGER(iwp) ::  nxrc            !<
+    INTEGER(iwp) ::  nxrf            !<
+    INTEGER(iwp) ::  nxr_on_file     !<
+    INTEGER(iwp) ::  nync            !<
+    INTEGER(iwp) ::  nynf            !<
+    INTEGER(iwp) ::  nyn_on_file     !<
+    INTEGER(iwp) ::  nysc            !<
+    INTEGER(iwp) ::  nysf            !<
+    INTEGER(iwp) ::  nys_on_file     !<
 
     LOGICAL, INTENT(OUT)  :: found
 
-    REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d   !< 
-
-    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !< 
-
-    REAL(wp), DIMENSION(0:0,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d2   !< 
+    REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d   !<
+
+    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
+
+    REAL(wp), DIMENSION(0:0,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d2   !<
 
 
@@ -12074 +12077 @@
           IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
              ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
-          ENDIF  
+          ENDIF
           IF ( k == 1 )  READ ( 13 )  tmp_2d
           rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =           &
                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-                        
+
        CASE ( 'rad_lw_in_xy_av' )
           IF ( .NOT. ALLOCATED( rad_lw_in_xy_av ) )  THEN
              ALLOCATE( rad_lw_in_xy_av(nysg:nyng,nxlg:nxrg) )
-          ENDIF  
+          ENDIF
           IF ( k == 1 )  READ ( 13 )  tmp_2d
           rad_lw_in_xy_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =      &
                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-                        
+
        CASE ( 'rad_lw_out_xy_av' )
           IF ( .NOT. ALLOCATED( rad_lw_out_xy_av ) )  THEN
              ALLOCATE( rad_lw_out_xy_av(nysg:nyng,nxlg:nxrg) )
-          ENDIF  
+          ENDIF
           IF ( k == 1 )  READ ( 13 )  tmp_2d
           rad_lw_out_xy_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =     &
                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-                        
+
        CASE ( 'rad_sw_in_xy_av' )
           IF ( .NOT. ALLOCATED( rad_sw_in_xy_av ) )  THEN
              ALLOCATE( rad_sw_in_xy_av(nysg:nyng,nxlg:nxrg) )
-          ENDIF  
+          ENDIF
           IF ( k == 1 )  READ ( 13 )  tmp_2d
           rad_sw_in_xy_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =      &
                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-                        
+
        CASE ( 'rad_sw_out_xy_av' )
           IF ( .NOT. ALLOCATED( rad_sw_out_xy_av ) )  THEN
              ALLOCATE( rad_sw_out_xy_av(nysg:nyng,nxlg:nxrg) )
-          ENDIF  
+          ENDIF
           IF ( k == 1 )  READ ( 13 )  tmp_2d
           rad_sw_out_xy_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =     &
                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-                        
+
        CASE ( 'rad_lw_in' )
           IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
@@ -12120 +12123 @@
                 ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
              ENDIF
-          ENDIF  
+          ENDIF
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
@@ -12144 +12147 @@
                 ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
              ENDIF
-          ENDIF  
+          ENDIF
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
@@ -12168 +12171 @@
                 ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
              ENDIF
-          ENDIF  
+          ENDIF
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
@@ -12192 +12195 @@
                 ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
              ENDIF
-          ENDIF  
+          ENDIF
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
@@ -12248 +12251 @@
                 ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
              ENDIF
-          ENDIF  
+          ENDIF
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
@@ -12272 +12275 @@
                 ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
              ENDIF
-          ENDIF  
+          ENDIF
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
@@ -12296 +12299 @@
                 ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
              ENDIF
-          ENDIF  
+          ENDIF
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
@@ -12320 +12323 @@
                 ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
              ENDIF
-          ENDIF  
+          ENDIF
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &