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

!> @file radiation_model_mod.f90
!------------------------------------------------------------------------------!
! This file is part of the PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software
! Foundation, either version 3 of the License, or (at your option) any later
! version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 2015-2018 Institute of Computer Science of the 
!                     Czech Academy of Sciences, Prague
! Copyright 2015-2018 Czech Technical University in Prague
! Copyright 1997-2019 Leibniz Universitaet Hannover
!------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
! 
! 
! Former revisions:
! -----------------
! $Id: radiation_model_mod.f90 3861 2019-04-04 06:27:41Z maronga $
! Bugfix: factor of 4.0 required instead of 3.0 in calculation of rad_lw_out_change_0
! 
! 3859 2019-04-03 20:30:31Z maronga
! Added some descriptions
! 
! 3847 2019-04-01 14:51:44Z suehring
! Implement check for dt_radiation (must be > 0)
! 
! 3846 2019-04-01 13:55:30Z suehring
! unused variable removed
! 
! 3814 2019-03-26 08:40:31Z pavelkrc
! Change zenith(0:0) and others to scalar.
! Code review.
! Rename exported nzu, nzp and related variables due to name conflict
! 
! 3771 2019-02-28 12:19:33Z raasch
! rrtmg preprocessor for directives moved/added, save attribute added to temporary
! pointers to avoid compiler warnings about outlived pointer targets,
! statement added to avoid compiler warning about unused variable
! 
! 3769 2019-02-28 10:16:49Z moh.hefny
! removed unused variables and subroutine radiation_radflux_gridbox
!
! 3767 2019-02-27 08:18:02Z raasch
! unused variable for file index removed from rrd-subroutines parameter list
! 
! 3760 2019-02-21 18:47:35Z moh.hefny
! Bugfix: initialized simulated_time before calculating solar position
! to enable restart option with reading in SVF from file(s).
!
! 3754 2019-02-19 17:02:26Z kanani
! (resler, pavelkrc)
! Bugfixes: add further required MRT factors to read/write_svf,
! fix for aggregating view factors to eliminate local noise in reflected 
! irradiance at mutually close surfaces (corners, presence of trees) in the 
! angular discretization scheme.
!
! 3752 2019-02-19 09:37:22Z resler
! added read/write number of MRT factors to the respective routines
!
! 3705 2019-01-29 19:56:39Z suehring
! Make variables that are sampled in virtual measurement module public 
! 
! 3704 2019-01-29 19:51:41Z suehring
! Some interface calls moved to module_interface + cleanup
! 
! 3667 2019-01-10 14:26:24Z schwenkel
! Modified check for rrtmg input files
! 
! 3655 2019-01-07 16:51:22Z knoop
! nopointer option removed
! 
! 3633 2018-12-17 16:17:57Z schwenkel
! Include check for rrtmg files
! 
! 3630 2018-12-17 11:04:17Z knoop
! - fix initialization of date and time after calling zenith
! - fix a bug in radiation_solar_pos
!
! 3616 2018-12-10 09:44:36Z Salim
! fix manipulation of time variables in radiation_presimulate_solar_pos
!
! 3608 2018-12-07 12:59:57Z suehring $
! Bugfix radiation output
! 
! 3607 2018-12-07 11:56:58Z suehring
! Output of radiation-related quantities migrated to radiation_model_mod.
! 
! 3589 2018-11-30 15:09:51Z suehring
! Remove erroneous UTF encoding
! 
! 3572 2018-11-28 11:40:28Z suehring
! Add filling the short- and longwave radiation flux arrays (e.g. diffuse, 
! direct, reflected, resedual) for all surfaces. This is required to surface
! outputs in suface_output_mod. (M. Salim)
! 
! 3571 2018-11-28 09:24:03Z moh.hefny
! Add an epsilon value to compare values in if statement to fix possible 
! precsion related errors in raytrace routines.
! 
! 3524 2018-11-14 13:36:44Z raasch
! missing cpp-directives added
! 
! 3495 2018-11-06 15:22:17Z kanani
! Resort control_parameters ONLY list,
! From branch radiation@3491 moh.hefny:
! bugfix in calculating the apparent solar positions by updating 
! the simulated time so that the actual time is correct.
! 
! 3464 2018-10-30 18:08:55Z kanani
! From branch resler@3462, pavelkrc:
! add MRT shaping function for human
! 
! 3449 2018-10-29 19:36:56Z suehring
! New RTM version 3.0: (Pavel Krc, Jaroslav Resler, ICS, Prague)
!   - Interaction of plant canopy with LW radiation
!   - Transpiration from resolved plant canopy dependent on radiation
!     called from RTM
! 
! 
! 3435 2018-10-26 18:25:44Z gronemeier
! - workaround: return unit=illegal in check_data_output for certain variables
!   when check called from init_masks
! - Use pointer in masked output to reduce code redundancies
! - Add terrain-following masked output
! 
! 3424 2018-10-25 07:29:10Z gronemeier
! bugfix: add rad_lw_in, rad_lw_out, rad_sw_out to radiation_check_data_output
! 
! 3378 2018-10-19 12:34:59Z kanani
! merge from radiation branch (r3362) into trunk
! (moh.hefny):
! - removed read/write_svf_on_init and read_dist_max_svf (not used anymore)
! - bugfix nzut > nzpt in calculating maxboxes
! 
! 3372 2018-10-18 14:03:19Z raasch
! bugfix: kind type of 2nd argument of mpi_win_allocate changed, misplaced
!         __parallel directive
! 
! 3351 2018-10-15 18:40:42Z suehring
! Do not overwrite values of spectral and broadband albedo during initialization
! if they are already initialized in the urban-surface model via ASCII input.
! 
! 3337 2018-10-12 15:17:09Z kanani
! - New RTM version 2.9: (Pavel Krc, Jaroslav Resler, ICS, Prague)
!   added calculation of the MRT inside the RTM module
!   MRT fluxes are consequently used in the new biometeorology module
!   for calculation of biological indices (MRT, PET)
!   Fixes of v. 2.5 and SVN trunk:
!    - proper initialization of rad_net_l
!    - make arrays nsurfs and surfstart TARGET to prevent some MPI problems
!    - initialization of arrays used in MPI one-sided operation as 1-D arrays
!      to prevent problems with some MPI/compiler combinations
!    - fix indexing of target displacement in subroutine request_itarget to 
!      consider nzub
!    - fix LAD dimmension range in PCB calculation
!    - check ierr in all MPI calls
!    - use proper per-gridbox sky and diffuse irradiance
!    - fix shading for reflected irradiance
!    - clear away the residuals of "atmospheric surfaces" implementation
!    - fix rounding bug in raytrace_2d introduced in SVN trunk
! - New RTM version 2.5: (Pavel Krc, Jaroslav Resler, ICS, Prague)
!   can use angular discretization for all SVF
!   (i.e. reflected and emitted radiation in addition to direct and diffuse),
!   allowing for much better scaling wih high resoltion and/or complex terrain
! - Unite array grow factors
! - Fix slightly shifted terrain height in raytrace_2d
! - Use more efficient MPI_Win_allocate for reverse gridsurf index
! - Fix random MPI RMA bugs on Intel compilers
! - Fix approx. double plant canopy sink values for reflected radiation
! - Fix mostly missing plant canopy sinks for direct radiation
! - Fix discretization errors for plant canopy sink in diffuse radiation
! - Fix rounding errors in raytrace_2d
!
! 3274 2018-09-24 15:42:55Z knoop
! Modularization of all bulk cloud physics code components
! 
! 3272 2018-09-24 10:16:32Z suehring
! - split direct and diffusion shortwave radiation using RRTMG rather than using
!   calc_diffusion_radiation, in case of RRTMG
! - removed the namelist variable split_diffusion_radiation. Now splitting depends
!   on the choise of radiation radiation scheme
! - removed calculating the rdiation flux for surfaces at the radiation scheme
!   in case of using RTM since it will be calculated anyway in the radiation 
!   interaction routine.
! - set SW radiation flux for surfaces to zero at night in case of no RTM is used
! - precalculate the unit vector yxdir of ray direction to avoid the temporarly 
!   array allocation during the subroutine call
! - fixed a bug in calculating the max number of boxes ray can cross in the domain
! 
! 3264 2018-09-20 13:54:11Z moh.hefny
! Bugfix in raytrace_2d calls
! 
! 3248 2018-09-14 09:42:06Z sward
! Minor formating changes
! 
! 3246 2018-09-13 15:14:50Z sward
! Added error handling for input namelist via parin_fail_message
! 
! 3241 2018-09-12 15:02:00Z raasch
! unused variables removed or commented
! 
! 3233 2018-09-07 13:21:24Z schwenkel
! Adapted for the use of cloud_droplets 
! 
! 3230 2018-09-05 09:29:05Z schwenkel
! Bugfix in radiation_constant_surf: changed (10.0 - emissivity_urb) to 
! (1.0 - emissivity_urb)
! 
! 3226 2018-08-31 12:27:09Z suehring
! Bugfixes in calculation of sky-view factors and canopy-sink factors.
! 
! 3186 2018-07-30 17:07:14Z suehring
! Remove print statement
! 
! 3180 2018-07-27 11:00:56Z suehring
! Revise concept for calculation of effective radiative temperature and mapping
! of radiative heating
! 
! 3175 2018-07-26 14:07:38Z suehring
! Bugfix for commit 3172
! 
! 3173 2018-07-26 12:55:23Z suehring
! Revise output of surface radiation quantities in case of overhanging 
! structures
! 
! 3172 2018-07-26 12:06:06Z suehring
! Bugfixes:
!  - temporal work-around for calculation of effective radiative surface 
!    temperature
!  - prevent positive solar radiation during nighttime
! 
! 3170 2018-07-25 15:19:37Z suehring
! Bugfix, map signle-column radiation forcing profiles on top of any topography
! 
! 3156 2018-07-19 16:30:54Z knoop
! Bugfix: replaced usage of the pt array with the surf%pt_surface array
!
! 3137 2018-07-17 06:44:21Z maronga
! String length for trace_names fixed
! 
! 3127 2018-07-15 08:01:25Z maronga
! A few pavement parameters updated.
! 
! 3123 2018-07-12 16:21:53Z suehring
! Correct working precision for INTEGER number
! 
! 3122 2018-07-11 21:46:41Z maronga
! Bugfix: maximum distance for raytracing was set to  -999 m by default,
! effectively switching off all surface reflections when max_raytracing_dist
! was not explicitly set in namelist
! 
! 3117 2018-07-11 09:59:11Z maronga
! Bugfix: water vapor was not transfered to RRTMG when bulk_cloud_model = .F.
! Bugfix: changed the calculation of RRTMG boundary conditions (Mohamed Salim)
! Bugfix: dry residual atmosphere is replaced by standard RRTMG atmosphere
! 
! 3116 2018-07-10 14:31:58Z suehring
! Output of long/shortwave radiation at surface
! 
! 3107 2018-07-06 15:55:51Z suehring
! Bugfix, missing index for dz
! 
! 3066 2018-06-12 08:55:55Z Giersch
! Error message revised 
! 
! 3065 2018-06-12 07:03:02Z Giersch
! dz was replaced by dz(1), error message concerning vertical stretching was 
! added  
! 
! 3049 2018-05-29 13:52:36Z Giersch
! Error messages revised
! 
! 3045 2018-05-28 07:55:41Z Giersch
! Error message revised
! 
! 3026 2018-05-22 10:30:53Z schwenkel
! Changed the name specific humidity to mixing ratio, since we are computing
! mixing ratios.
! 
! 3016 2018-05-09 10:53:37Z Giersch
! Revised structure of reading svf data according to PALM coding standard: 
! svf_code_field/len and fsvf removed, error messages PA0493 and PA0494 added,
! allocation status of output arrays checked. 
! 
! 3014 2018-05-09 08:42:38Z maronga
! Introduced plant canopy height similar to urban canopy height to limit
! the memory requirement to allocate lad.
! Deactivated automatic setting of minimum raytracing distance.
! 
! 3004 2018-04-27 12:33:25Z Giersch
! Further allocation checks implemented (averaged data will be assigned to fill
! values if no allocation happened so far) 
! 
! 2995 2018-04-19 12:13:16Z Giersch
! IF-statement in radiation_init removed so that the calculation of radiative 
! fluxes at model start is done in any case, bugfix in 
! radiation_presimulate_solar_pos (end_time is the sum of end_time and the 
! spinup_time specified in the p3d_file ), list of variables/fields that have
! to be written out or read in case of restarts has been extended
! 
! 2977 2018-04-17 10:27:57Z kanani
! Implement changes from branch radiation (r2948-2971) with minor modifications,
! plus some formatting.
! (moh.hefny):
! - replaced plant_canopy by npcbl to check tree existence to avoid weird 
!   allocation of related arrays (after domain decomposition some domains 
!   contains no trees although plant_canopy (global parameter) is still TRUE).
! - added a namelist parameter to force RTM settings
! - enabled the option to switch radiation reflections off
! - renamed surf_reflections to surface_reflections 
! - removed average_radiation flag from the namelist (now it is implicitly set 
!   in init_3d_model according to RTM)
! - edited read and write sky view factors and CSF routines to account for
!   the sub-domains which may not contain any of them
! 
! 2967 2018-04-13 11:22:08Z raasch
! bugfix: missing parallel cpp-directives added
! 
! 2964 2018-04-12 16:04:03Z Giersch
! Error message PA0491 has been introduced which could be previously found in 
! check_open. The variable numprocs_previous_run is only known in case of 
! initializing_actions == read_restart_data
! 
! 2963 2018-04-12 14:47:44Z suehring
! - Introduce index for vegetation/wall, pavement/green-wall and water/window 
!   surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
! - Minor bugfix in initialization of albedo for window surfaces
! 
! 2944 2018-04-03 16:20:18Z suehring
! Fixed bad commit
! 
! 2943 2018-04-03 16:17:10Z suehring
! No read of nsurfl from SVF file since it is calculated in
! radiation_interaction_init,
! allocation of arrays in radiation_read_svf only if not yet allocated,
! update of 2920 revision comment.
! 
! 2932 2018-03-26 09:39:22Z maronga
! renamed radiation_par to radiation_parameters
! 
! 2930 2018-03-23 16:30:46Z suehring
! Remove default surfaces from radiation model, does not make much sense to 
! apply radiation model without energy-balance solvers; Further, add check for 
! this. 
! 
! 2920 2018-03-22 11:22:01Z kanani
! - Bugfix: Initialize pcbl array (=-1)
! RTM version 2.0 (Jaroslav Resler, Pavel Krc, Mohamed Salim):
! - new major version of radiation interactions
! - substantially enhanced performance and scalability
! - processing of direct and diffuse solar radiation separated from reflected 
!   radiation, removed virtual surfaces
! - new type of sky discretization by azimuth and elevation angles
! - diffuse radiation processed cumulatively using sky view factor
! - used precalculated apparent solar positions for direct irradiance
! - added new 2D raytracing process for processing whole vertical column at once
!   to increase memory efficiency and decrease number of MPI RMA operations
! - enabled limiting the number of view factors between surfaces by the distance
!   and value
! - fixing issues induced by transferring radiation interactions from 
!   urban_surface_mod to radiation_mod
! - bugfixes and other minor enhancements
! 
! 2906 2018-03-19 08:56:40Z Giersch
! NAMELIST paramter read/write_svf_on_init have been removed, functions 
! check_open and close_file are used now for opening/closing files related to
! svf data, adjusted unit number and error numbers 
! 
! 2894 2018-03-15 09:17:58Z Giersch
! Calculations of the index range of the subdomain on file which overlaps with
! the current subdomain are already done in read_restart_data_mod
! radiation_read_restart_data was renamed to radiation_rrd_local and 
! radiation_last_actions was renamed to radiation_wrd_local, variable named 
! found has been introduced for checking if restart data was found, reading 
! of restart strings has been moved completely to read_restart_data_mod, 
! radiation_rrd_local is already inside the overlap loop programmed in 
! read_restart_data_mod, the marker *** end rad *** is not necessary anymore,
! strings and their respective lengths are written out and read now in case of 
! restart runs to get rid of prescribed character lengths (Giersch)
!
! 2809 2018-02-15 09:55:58Z suehring
! Bugfix for gfortran: Replace the function C_SIZEOF with STORAGE_SIZE
! 
! 2753 2018-01-16 14:16:49Z suehring
! Tile approach for spectral albedo implemented. 
! 
! 2746 2018-01-15 12:06:04Z suehring
! Move flag plant canopy to modules
! 
! 2724 2018-01-05 12:12:38Z maronga
! Set default of average_radiation to .FALSE.
! 
! 2723 2018-01-05 09:27:03Z maronga
! Bugfix in calculation of rad_lw_out (clear-sky). First grid level was used
! instead of the surface value
! 
! 2718 2018-01-02 08:49:38Z maronga
! Corrected "Former revisions" section
! 
! 2707 2017-12-18 18:34:46Z suehring
! Changes from last commit documented
! 
! 2706 2017-12-18 18:33:49Z suehring
! Bugfix, in average radiation case calculate exner function before using it.
!
! 2701 2017-12-15 15:40:50Z suehring
! Changes from last commit documented
! 
! 2698 2017-12-14 18:46:24Z suehring
! Bugfix in get_topography_top_index
!
! 2696 2017-12-14 17:12:51Z kanani
! - Change in file header (GPL part)
! - Improved reading/writing of SVF from/to file (BM)
! - Bugfixes concerning RRTMG as well as average_radiation options (M. Salim)
! - Revised initialization of surface albedo and some minor bugfixes (MS)
! - Update net radiation after running radiation interaction routine (MS)
! - Revisions from M Salim included
! - Adjustment to topography and surface structure (MS)
! - Initialization of albedo and surface emissivity via input file (MS)
! - albedo_pars extended (MS)
! 
! 2604 2017-11-06 13:29:00Z schwenkel
! bugfix for calculation of effective radius using morrison microphysics
! 
! 2601 2017-11-02 16:22:46Z scharf
! added emissivity to namelist
! 
! 2575 2017-10-24 09:57:58Z maronga
! Bugfix: calculation of shortwave and longwave albedos for RRTMG swapped
! 
! 2547 2017-10-16 12:41:56Z schwenkel
! extended by cloud_droplets option, minor bugfix and correct calculation of
! cloud droplet number concentration
! 
! 2544 2017-10-13 18:09:32Z maronga
! Moved date and time quantitis to separate module date_and_time_mod
! 
! 2512 2017-10-04 08:26:59Z raasch
! upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny
! no output of ghost layer data
! 
! 2504 2017-09-27 10:36:13Z maronga
! Updates pavement types and albedo parameters
! 
! 2328 2017-08-03 12:34:22Z maronga
! Emissivity can now be set individually for each pixel.
! Albedo type can be inferred from land surface model.
! Added default albedo type for bare soil
! 
! 2318 2017-07-20 17:27:44Z suehring
! Get topography top index via Function call 
! 
! 2317 2017-07-20 17:27:19Z suehring
! Improved syntax layout
! 
! 2298 2017-06-29 09:28:18Z raasch
! type of write_binary changed from CHARACTER to LOGICAL
! 
! 2296 2017-06-28 07:53:56Z maronga
! Added output of rad_sw_out for radiation_scheme = 'constant'
! 
! 2270 2017-06-09 12:18:47Z maronga
! Numbering changed (2 timeseries removed)
! 
! 2249 2017-06-06 13:58:01Z sward
! Allow for RRTMG runs without humidity/cloud physics
! 
! 2248 2017-06-06 13:52:54Z sward
! Error no changed
! 
! 2233 2017-05-30 18:08:54Z suehring
! 
! 2232 2017-05-30 17:47:52Z suehring
! Adjustments to new topography concept
! Bugfix in read restart
! 
! 2200 2017-04-11 11:37:51Z suehring
! Bugfix in call of exchange_horiz_2d and read restart data
! 
! 2163 2017-03-01 13:23:15Z schwenkel
! Bugfix in radiation_check_data_output
! 
! 2157 2017-02-22 15:10:35Z suehring
! Bugfix in read_restart data
! 
! 2011 2016-09-19 17:29:57Z kanani
! Removed CALL of auxiliary SUBROUTINE get_usm_info,
! flag urban_surface is now defined in module control_parameters.
! 
! 2007 2016-08-24 15:47:17Z kanani
! Added calculation of solar directional vector for new urban surface 
! model,
! accounted for urban_surface model in radiation_check_parameters,
! correction of comments for zenith angle.
! 
! 2000 2016-08-20 18:09:15Z knoop
! Forced header and separation lines into 80 columns
! 
! 1976 2016-07-27 13:28:04Z maronga
! Output of 2D/3D/masked data is now directly done within this module. The
! radiation schemes have been simplified for better usability so that
! rad_lw_in, rad_lw_out, rad_sw_in, and rad_sw_out are available independent of
! the radiation code used.
! 
! 1856 2016-04-13 12:56:17Z maronga
! Bugfix: allocation of rad_lw_out for radiation_scheme = 'clear-sky'
! 
! 1853 2016-04-11 09:00:35Z maronga
! Added routine for radiation_scheme = constant.
!  
! 1849 2016-04-08 11:33:18Z hoffmann 
! Adapted for modularization of microphysics
!
! 1826 2016-04-07 12:01:39Z maronga
! Further modularization.
! 
! 1788 2016-03-10 11:01:04Z maronga
! Added new albedo class for pavements / roads.
!
! 1783 2016-03-06 18:36:17Z raasch
! palm-netcdf-module removed in order to avoid a circular module dependency,
! netcdf-variables moved to netcdf-module, new routine netcdf_handle_error_rad
! added
!
! 1757 2016-02-22 15:49:32Z maronga
! Added parameter unscheduled_radiation_calls. Bugfix: interpolation of sounding
! profiles for pressure and temperature above the LES domain.
! 
! 1709 2015-11-04 14:47:01Z maronga
! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
! corrections
! 
! 1701 2015-11-02 07:43:04Z maronga
! Bugfixes: wrong index for output of timeseries, setting of nz_snd_end
! 
! 1691 2015-10-26 16:17:44Z maronga
! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
! Added output of radiative heating rates.
! 
! 1682 2015-10-07 23:56:08Z knoop
! Code annotations made doxygen readable 
! 
! 1606 2015-06-29 10:43:37Z maronga
! Added preprocessor directive __netcdf to allow for compiling without netCDF.
! Note, however, that RRTMG cannot be used without netCDF.
! 
! 1590 2015-05-08 13:56:27Z maronga
! Bugfix: definition of character strings requires same length for all elements
! 
! 1587 2015-05-04 14:19:01Z maronga
! Added albedo class for snow
! 
! 1585 2015-04-30 07:05:52Z maronga
! Added support for RRTMG
! 
! 1571 2015-03-12 16:12:49Z maronga
! Added missing KIND attribute. Removed upper-case variable names
! 
! 1551 2015-03-03 14:18:16Z maronga
! Added support for data output. Various variables have been renamed. Added
! interface for different radiation schemes (currently: clear-sky, constant, and
! RRTM (not yet implemented).
! 
! 1496 2014-12-02 17:25:50Z maronga
! Initial revision
! 
!
! Description:
! ------------
!> Radiation models and interfaces
!> @todo Replace dz(1) appropriatly to account for grid stretching
!> @todo move variable definitions used in radiation_init only to the subroutine
!>       as they are no longer required after initialization.
!> @todo Output of full column vertical profiles used in RRTMG
!> @todo Output of other rrtm arrays (such as volume mixing ratios)
!> @todo Check for mis-used NINT() calls in raytrace_2d
!>       RESULT: Original was correct (carefully verified formula), the change
!>               to INT broke raytracing      -- P. Krc
!> @todo Optimize radiation_tendency routines
!>
!> @note Many variables have a leading dummy dimension (0:0) in order to
!>       match the assume-size shape expected by the RRTMG model.
!------------------------------------------------------------------------------!
 MODULE radiation_model_mod
 
    USE arrays_3d,                                                             &
        ONLY:  dzw, hyp, nc, pt, p, q, ql, u, v, w, zu, zw, exner, d_exner

    USE basic_constants_and_equations_mod,                                     &
        ONLY:  c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, sigma_sb, &
               barometric_formula

    USE calc_mean_profile_mod,                                                 &
        ONLY:  calc_mean_profile

    USE control_parameters,                                                    &
        ONLY:  cloud_droplets, coupling_char, dz, dt_spinup, end_time,         &
               humidity,                                                       &
               initializing_actions, io_blocks, io_group,                      &
               land_surface, large_scale_forcing,                              &
               latitude, longitude, lsf_surf,                                  &
               message_string, plant_canopy, pt_surface,                       &
               rho_surface, simulated_time, spinup_time, surface_pressure,     &
               read_svf, write_svf,                                            &
               time_since_reference_point, urban_surface, varnamelength

    USE cpulog,                                                                &
        ONLY:  cpu_log, log_point, log_point_s

    USE grid_variables,                                                        &
         ONLY:  ddx, ddy, dx, dy 

    USE date_and_time_mod,                                                     &
        ONLY:  calc_date_and_time, d_hours_day, d_seconds_hour, d_seconds_year,&
               day_of_year, d_seconds_year, day_of_month, day_of_year_init,    &
               init_date_and_time, month_of_year, time_utc_init, time_utc

    USE indices,                                                               &
        ONLY:  nnx, nny, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,   &
               nzb, nzt

    USE, INTRINSIC :: iso_c_binding

    USE kinds

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

#if defined ( __netcdf )
    USE NETCDF
#endif

    USE netcdf_data_input_mod,                                                 &
        ONLY:  albedo_type_f, albedo_pars_f, building_type_f, pavement_type_f, &
               vegetation_type_f, water_type_f

    USE plant_canopy_model_mod,                                                &
        ONLY:  lad_s, pc_heating_rate, pc_transpiration_rate, pc_latent_rate,  &
               plant_canopy_transpiration, pcm_calc_transpiration_rate

    USE pegrid

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

    USE parrrtm,                                                               &
        ONLY:  nbndlw

    USE rrtmg_lw_init,                                                         &
        ONLY:  rrtmg_lw_ini

    USE rrtmg_sw_init,                                                         &
        ONLY:  rrtmg_sw_ini

    USE rrtmg_lw_rad,                                                          &
        ONLY:  rrtmg_lw

    USE rrtmg_sw_rad,                                                          &
        ONLY:  rrtmg_sw
#endif
    USE statistics,                                                            &
        ONLY:  hom

    USE surface_mod,                                                           &
        ONLY:  get_topography_top_index, get_topography_top_index_ji,          &
               ind_pav_green, ind_veg_wall, ind_wat_win,                       &
               surf_lsm_h, surf_lsm_v, surf_type, surf_usm_h, surf_usm_v,      &
               vertical_surfaces_exist

    IMPLICIT NONE

    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'

!
!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
    CHARACTER(37), DIMENSION(0:33), PARAMETER :: albedo_type_name = (/      &
                                   'user defined                         ', & !  0 
                                   'ocean                                ', & !  1
                                   'mixed farming, tall grassland        ', & !  2
                                   'tall/medium grassland                ', & !  3 
                                   'evergreen shrubland                  ', & !  4
                                   'short grassland/meadow/shrubland     ', & !  5
                                   'evergreen needleleaf forest          ', & !  6
                                   'mixed deciduous evergreen forest     ', & !  7
                                   'deciduous forest                     ', & !  8
                                   'tropical evergreen broadleaved forest', & !  9
                                   'medium/tall grassland/woodland       ', & ! 10
                                   'desert, sandy                        ', & ! 11 
                                   'desert, rocky                        ', & ! 12 
                                   'tundra                               ', & ! 13
                                   'land ice                             ', & ! 14 
                                   'sea ice                              ', & ! 15 
                                   'snow                                 ', & ! 16
                                   'bare soil                            ', & ! 17
                                   'asphalt/concrete mix                 ', & ! 18
                                   'asphalt (asphalt concrete)           ', & ! 19
                                   'concrete (Portland concrete)         ', & ! 20
                                   'sett                                 ', & ! 21
                                   'paving stones                        ', & ! 22
                                   'cobblestone                          ', & ! 23
                                   'metal                                ', & ! 24
                                   'wood                                 ', & ! 25
                                   'gravel                               ', & ! 26
                                   'fine gravel                          ', & ! 27
                                   'pebblestone                          ', & ! 28
                                   'woodchips                            ', & ! 29
                                   'tartan (sports)                      ', & ! 30
                                   'artifical turf (sports)              ', & ! 31
                                   'clay (sports)                        ', & ! 32
                                   'building (dummy)                     '  & ! 33
                                                         /)

    INTEGER(iwp) :: albedo_type  = 9999999_iwp, &     !< Albedo surface type
                    dots_rad     = 0_iwp              !< starting index for timeseries output 

    LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
                constant_albedo = .FALSE.,            & !< flag parameter indicating whether the albedo may change depending on zenith
                force_radiation_call = .FALSE.,       & !< flag parameter for unscheduled radiation calls
                lw_radiation = .TRUE.,                & !< flag parameter indicating whether longwave radiation shall be calculated
                radiation = .FALSE.,                  & !< flag parameter indicating whether the radiation model is used
                sun_up    = .TRUE.,                   & !< flag parameter indicating whether the sun is up or down
                sw_radiation = .TRUE.,                & !< flag parameter indicating whether shortwave radiation shall be calculated
                sun_direction = .FALSE.,              & !< flag parameter indicating whether solar direction shall be calculated
                average_radiation = .FALSE.,          & !< flag to set the calculation of radiation averaging for the domain
                radiation_interactions = .FALSE.,     & !< flag to activiate RTM (TRUE only if vertical urban/land surface and trees exist)
                surface_reflections = .TRUE.,         & !< flag to switch the calculation of radiation interaction between surfaces.
                                                        !< When it switched off, only the effect of buildings and trees shadow
                                                        !< will be considered. However fewer SVFs are expected.
                radiation_interactions_on = .TRUE.      !< namelist flag to force RTM activiation regardless to vertical urban/land surface and trees

    REAL(wp) :: albedo = 9999999.9_wp,           & !< NAMELIST alpha
                albedo_lw_dif = 9999999.9_wp,    & !< NAMELIST aldif
                albedo_lw_dir = 9999999.9_wp,    & !< NAMELIST aldir
                albedo_sw_dif = 9999999.9_wp,    & !< NAMELIST asdif
                albedo_sw_dir = 9999999.9_wp,    & !< NAMELIST asdir
                decl_1,                          & !< declination coef. 1
                decl_2,                          & !< declination coef. 2
                decl_3,                          & !< declination coef. 3
                dt_radiation = 0.0_wp,           & !< radiation model timestep
                emissivity = 9999999.9_wp,       & !< NAMELIST surface emissivity
                lon = 0.0_wp,                    & !< longitude in radians
                lat = 0.0_wp,                    & !< latitude in radians
                net_radiation = 0.0_wp,          & !< net radiation at surface
                skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
                sky_trans,                       & !< sky transmissivity
                time_radiation = 0.0_wp            !< time since last call of radiation code


    REAL(wp) ::  cos_zenith        !< cosine of solar zenith angle, also z-coordinate of solar unit vector
    REAL(wp) ::  sun_dir_lat       !< y-coordinate of solar unit vector
    REAL(wp) ::  sun_dir_lon       !< x-coordinate of solar unit vector

    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_net_av       !< average of net radiation (rad_net) at surface
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_in_xy_av  !< average of incoming longwave radiation at surface
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_out_xy_av !< average of outgoing longwave radiation at surface
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_in_xy_av  !< average of incoming shortwave radiation at surface
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_out_xy_av !< average of outgoing shortwave radiation at surface

    REAL(wp), PARAMETER :: emissivity_atm_clsky = 0.8_wp       !< emissivity of the clear-sky atmosphere
!
!-- Land surface albedos for solar zenith angle of 60degree after Briegleb (1992)     
!-- (shortwave, longwave, broadband):   sw,      lw,      bb,
    REAL(wp), DIMENSION(0:2,1:33), PARAMETER :: albedo_pars = RESHAPE( (/& 
                                   0.06_wp, 0.06_wp, 0.06_wp,            & !  1
                                   0.09_wp, 0.28_wp, 0.19_wp,            & !  2
                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  3
                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  4
                                   0.14_wp, 0.34_wp, 0.25_wp,            & !  5
                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  6
                                   0.06_wp, 0.27_wp, 0.17_wp,            & !  7
                                   0.06_wp, 0.31_wp, 0.19_wp,            & !  8
                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  9
                                   0.06_wp, 0.28_wp, 0.18_wp,            & ! 10
                                   0.35_wp, 0.51_wp, 0.43_wp,            & ! 11
                                   0.24_wp, 0.40_wp, 0.32_wp,            & ! 12
                                   0.10_wp, 0.27_wp, 0.19_wp,            & ! 13
                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 14
                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 15
                                   0.95_wp, 0.70_wp, 0.82_wp,            & ! 16
                                   0.08_wp, 0.08_wp, 0.08_wp,            & ! 17
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 18
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 19
                                   0.30_wp, 0.30_wp, 0.30_wp,            & ! 20
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 21
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 22
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 23
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 24
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 25
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 26
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 27
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 28
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 29
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 30
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 31
                                   0.17_wp, 0.17_wp, 0.17_wp,            & ! 32
                                   0.17_wp, 0.17_wp, 0.17_wp             & ! 33
                                 /), (/ 3, 33 /) )

    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
                        rad_lw_cs_hr,                  & !< longwave clear sky radiation heating rate (K/s)
                        rad_lw_cs_hr_av,               & !< average of rad_lw_cs_hr
                        rad_lw_hr,                     & !< longwave radiation heating rate (K/s)
                        rad_lw_hr_av,                  & !< average of rad_sw_hr
                        rad_lw_in,                     & !< incoming longwave radiation (W/m2)
                        rad_lw_in_av,                  & !< average of rad_lw_in
                        rad_lw_out,                    & !< outgoing longwave radiation (W/m2)
                        rad_lw_out_av,                 & !< average of rad_lw_out
                        rad_sw_cs_hr,                  & !< shortwave clear sky radiation heating rate (K/s)
                        rad_sw_cs_hr_av,               & !< average of rad_sw_cs_hr
                        rad_sw_hr,                     & !< shortwave radiation heating rate (K/s)
                        rad_sw_hr_av,                  & !< average of rad_sw_hr
                        rad_sw_in,                     & !< incoming shortwave radiation (W/m2)
                        rad_sw_in_av,                  & !< average of rad_sw_in
                        rad_sw_out,                    & !< outgoing shortwave radiation (W/m2)
                        rad_sw_out_av                    !< average of rad_sw_out


!
!-- Variables and parameters used in RRTMG only
#if defined ( __rrtmg )
    CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)


!
!-- Flag parameters to be passed to RRTMG (should not be changed until ice phase in clouds is allowed)
    INTEGER(iwp), PARAMETER :: rrtm_idrv     = 1, & !< flag for longwave upward flux calculation option (0,1)
                               rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
                               rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
                               rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
                               rrtm_inflgsw  = 2, & !< flag for sw cloud optical properties (0,1,2)
                               rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
                               rrtm_liqflgsw = 1    !< flag for sw liquid droplet specifications

!
!-- The following variables should be only changed with care, as this will 
!-- require further setting of some variables, which is currently not
!-- implemented (aerosols, ice phase).
    INTEGER(iwp) :: nzt_rad,           & !< upper vertical limit for radiation calculations
                    rrtm_icld = 0,     & !< cloud flag (0: clear sky column, 1: cloudy column)
                    rrtm_iaer = 0        !< aerosol option flag (0: no aerosol layers, for lw only: 6 (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented)

    INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling

    LOGICAL :: snd_exists = .FALSE.      !< flag parameter to check whether a user-defined input files exists
    LOGICAL :: sw_exists = .FALSE.       !< flag parameter to check whether that required rrtmg sw file exists
    LOGICAL :: lw_exists = .FALSE.       !< flag parameter to check whether that required rrtmg lw file exists


    REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor

    REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !< hypostatic pressure from sounding data (hPa)
                                           rrtm_tsfc,   & !< dummy array for storing surface temperature
                                           t_snd          !< actual temperature from sounding data (hPa)

    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rrtm_ccl4vmr,   & !< CCL4 volume mixing ratio (g/mol)
                                             rrtm_cfc11vmr,  & !< CFC11 volume mixing ratio (g/mol)
                                             rrtm_cfc12vmr,  & !< CFC12 volume mixing ratio (g/mol)
                                             rrtm_cfc22vmr,  & !< CFC22 volume mixing ratio (g/mol)
                                             rrtm_ch4vmr,    & !< CH4 volume mixing ratio
                                             rrtm_cicewp,    & !< in-cloud ice water path (g/m2)
                                             rrtm_cldfr,     & !< cloud fraction (0,1)
                                             rrtm_cliqwp,    & !< in-cloud liquid water path (g/m2)
                                             rrtm_co2vmr,    & !< CO2 volume mixing ratio (g/mol)
                                             rrtm_emis,      & !< surface emissivity (0-1)  
                                             rrtm_h2ovmr,    & !< H2O volume mixing ratio
                                             rrtm_n2ovmr,    & !< N2O volume mixing ratio
                                             rrtm_o2vmr,     & !< O2 volume mixing ratio
                                             rrtm_o3vmr,     & !< O3 volume mixing ratio
                                             rrtm_play,      & !< pressure layers (hPa, zu-grid)
                                             rrtm_plev,      & !< pressure layers (hPa, zw-grid)
                                             rrtm_reice,     & !< cloud ice effective radius (microns)
                                             rrtm_reliq,     & !< cloud water drop effective radius (microns)
                                             rrtm_tlay,      & !< actual temperature (K, zu-grid)
                                             rrtm_tlev,      & !< actual temperature (K, zw-grid)
                                             rrtm_lwdflx,    & !< RRTM output of incoming longwave radiation flux (W/m2)
                                             rrtm_lwdflxc,   & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2) 
                                             rrtm_lwuflx,    & !< RRTM output of outgoing longwave radiation flux (W/m2)
                                             rrtm_lwuflxc,   & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
                                             rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
                                             rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
                                             rrtm_lwhr,      & !< RRTM output of longwave radiation heating rate (K/d)
                                             rrtm_lwhrc,     & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
                                             rrtm_swdflx,    & !< RRTM output of incoming shortwave radiation flux (W/m2)
                                             rrtm_swdflxc,   & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2) 
                                             rrtm_swuflx,    & !< RRTM output of outgoing shortwave radiation flux (W/m2)
                                             rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
                                             rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
                                             rrtm_swhrc,     & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d) 
                                             rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m2) 
                                             rrtm_difdflux     !< RRTM output of incoming diffuse shortwave (W/m2) 

    REAL(wp), DIMENSION(1) ::                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
                                             rrtm_aldir,     & !< surface albedo for longwave direct radiation
                                             rrtm_asdif,     & !< surface albedo for shortwave diffuse radiation
                                             rrtm_asdir        !< surface albedo for shortwave direct radiation

!
!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   & !< incoming clear sky longwave radiation (W/m2) (not used)
                                                rad_lw_cs_out,  & !< outgoing clear sky longwave radiation (W/m2) (not used)
                                                rad_sw_cs_in,   & !< incoming clear sky shortwave radiation (W/m2) (not used)
                                                rad_sw_cs_out,  & !< outgoing clear sky shortwave radiation (W/m2) (not used)
                                                rrtm_lw_tauaer, & !< lw aerosol optical depth
                                                rrtm_lw_taucld, & !< lw in-cloud optical depth
                                                rrtm_sw_taucld, & !< sw in-cloud optical depth
                                                rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
                                                rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
                                                rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
                                                rrtm_sw_tauaer, & !< sw aerosol optical depth
                                                rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
                                                rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
                                                rrtm_sw_ecaer     !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)

#endif
!
!-- Parameters of urban and land surface models
    INTEGER(iwp)                                   ::  nz_urban                           !< number of layers of urban surface (will be calculated)
    INTEGER(iwp)                                   ::  nz_plant                           !< number of layers of plant canopy (will be calculated)
    INTEGER(iwp)                                   ::  nz_urban_b                         !< bottom layer of urban surface (will be calculated)
    INTEGER(iwp)                                   ::  nz_urban_t                         !< top layer of urban surface (will be calculated)
    INTEGER(iwp)                                   ::  nz_plant_t                         !< top layer of plant canopy (will be calculated)
!-- parameters of urban and land surface models
    INTEGER(iwp), PARAMETER                        ::  nzut_free = 3                      !< number of free layers above top of of topography
    INTEGER(iwp), PARAMETER                        ::  ndsvf = 2                          !< number of dimensions of real values in SVF
    INTEGER(iwp), PARAMETER                        ::  idsvf = 2                          !< number of dimensions of integer values in SVF
    INTEGER(iwp), PARAMETER                        ::  ndcsf = 1                          !< number of dimensions of real values in CSF
    INTEGER(iwp), PARAMETER                        ::  idcsf = 2                          !< number of dimensions of integer values in CSF
    INTEGER(iwp), PARAMETER                        ::  kdcsf = 4                          !< number of dimensions of integer values in CSF calculation array
    INTEGER(iwp), PARAMETER                        ::  id = 1                             !< position of d-index in surfl and surf
    INTEGER(iwp), PARAMETER                        ::  iz = 2                             !< position of k-index in surfl and surf
    INTEGER(iwp), PARAMETER                        ::  iy = 3                             !< position of j-index in surfl and surf
    INTEGER(iwp), PARAMETER                        ::  ix = 4                             !< position of i-index in surfl and surf
    INTEGER(iwp), PARAMETER                        ::  im = 5                             !< position of surface m-index in surfl and surf
    INTEGER(iwp), PARAMETER                        ::  nidx_surf = 5                      !< number of indices in surfl and surf

    INTEGER(iwp), PARAMETER                        ::  nsurf_type = 10                    !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1

    INTEGER(iwp), PARAMETER                        ::  iup_u    = 0                       !< 0 - index of urban upward surface (ground or roof)
    INTEGER(iwp), PARAMETER                        ::  idown_u  = 1                       !< 1 - index of urban downward surface (overhanging)
    INTEGER(iwp), PARAMETER                        ::  inorth_u = 2                       !< 2 - index of urban northward facing wall
    INTEGER(iwp), PARAMETER                        ::  isouth_u = 3                       !< 3 - index of urban southward facing wall
    INTEGER(iwp), PARAMETER                        ::  ieast_u  = 4                       !< 4 - index of urban eastward facing wall
    INTEGER(iwp), PARAMETER                        ::  iwest_u  = 5                       !< 5 - index of urban westward facing wall

    INTEGER(iwp), PARAMETER                        ::  iup_l    = 6                       !< 6 - index of land upward surface (ground or roof)
    INTEGER(iwp), PARAMETER                        ::  inorth_l = 7                       !< 7 - index of land northward facing wall
    INTEGER(iwp), PARAMETER                        ::  isouth_l = 8                       !< 8 - index of land southward facing wall
    INTEGER(iwp), PARAMETER                        ::  ieast_l  = 9                       !< 9 - index of land eastward facing wall
    INTEGER(iwp), PARAMETER                        ::  iwest_l  = 10                      !< 10- index of land westward facing wall

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


!-- indices and sizes of urban and land surface models
    INTEGER(iwp)                                   ::  startland        !< start index of block of land and roof surfaces
    INTEGER(iwp)                                   ::  endland          !< end index of block of land and roof surfaces
    INTEGER(iwp)                                   ::  nlands           !< number of land and roof surfaces in local processor
    INTEGER(iwp)                                   ::  startwall        !< start index of block of wall surfaces
    INTEGER(iwp)                                   ::  endwall          !< end index of block of wall surfaces
    INTEGER(iwp)                                   ::  nwalls           !< number of wall surfaces in local processor

!-- indices needed for RTM netcdf output subroutines
    INTEGER(iwp), PARAMETER                        :: nd = 5
    CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_u = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_l = (/ iup_l, isouth_l, inorth_l, iwest_l, ieast_l /)
    INTEGER(iwp), DIMENSION(0:nd-1)                :: dirstart
    INTEGER(iwp), DIMENSION(0:nd-1)                :: dirend

!-- indices and sizes of urban and land surface models
    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surfl            !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x, m]
    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_linear     !< dtto (linearly allocated array)
    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surf             !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x, m]
    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surf_linear      !< dtto (linearly allocated array)
    INTEGER(iwp)                                   ::  nsurfl           !< number of all surfaces in local processor
    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  nsurfs           !< array of number of all surfaces in individual processors
    INTEGER(iwp)                                   ::  nsurf            !< global number of surfaces in index array of surfaces (nsurf = proc nsurfs)
    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfstart        !< starts of blocks of surfaces for individual processors in array surf (indexed from 1)
                                                                        !< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1]

!-- block variables needed for calculation of the plant canopy model inside the urban surface model
    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pct              !< top layer of the plant canopy
    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pch              !< heights of the plant canopy
    INTEGER(iwp)                                   ::  npcbl = 0        !< number of the plant canopy gridboxes in local processor
    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pcbl             !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i]
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw          !< array of absorbed sw radiation for local plant canopy box
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir       !< array of absorbed direct sw radiation for local plant canopy box
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif       !< array of absorbed diffusion sw radiation for local plant canopy box
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw          !< array of absorbed lw radiation for local plant canopy box

!-- configuration parameters (they can be setup in PALM config)
    LOGICAL                                        ::  raytrace_mpi_rma = .TRUE.          !< use MPI RMA to access LAD and gridsurf from remote processes during raytracing
    LOGICAL                                        ::  rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for
                                                                                          !< reflected radiation (as opposed to all mutually visible pairs)
    LOGICAL                                        ::  plant_lw_interact = .TRUE.         !< whether plant canopy interacts with LW radiation (in addition to SW)
    INTEGER(iwp)                                   ::  mrt_nlevels = 0                    !< number of vertical boxes above surface for which to calculate MRT
    LOGICAL                                        ::  mrt_skip_roof = .TRUE.             !< do not calculate MRT above roof surfaces
    LOGICAL                                        ::  mrt_include_sw = .TRUE.            !< should MRT calculation include SW radiation as well?
    LOGICAL                                        ::  mrt_geom_human = .TRUE.            !< MRT direction weights simulate human body instead of a sphere
    INTEGER(iwp)                                   ::  nrefsteps = 3                      !< number of reflection steps to perform
    REAL(wp), PARAMETER                            ::  ext_coef = 0.6_wp                  !< extinction coefficient (a.k.a. alpha)
    INTEGER(iwp), PARAMETER                        ::  rad_version_len = 10               !< length of identification string of rad version
    CHARACTER(rad_version_len), PARAMETER          ::  rad_version = 'RAD v. 3.0'         !< identification of version of binary svf and restart files
    INTEGER(iwp)                                   ::  raytrace_discrete_elevs = 40       !< number of discretization steps for elevation (nadir to zenith)
    INTEGER(iwp)                                   ::  raytrace_discrete_azims = 80       !< number of discretization steps for azimuth (out of 360 degrees)
    REAL(wp)                                       ::  max_raytracing_dist = -999.0_wp    !< maximum distance for raytracing (in metres)
    REAL(wp)                                       ::  min_irrf_value = 1e-6_wp           !< minimum potential irradiance factor value for raytracing
    REAL(wp), DIMENSION(1:30)                      ::  svfnorm_report_thresh = 1e21_wp    !< thresholds of SVF normalization values to report
    INTEGER(iwp)                                   ::  svfnorm_report_num                 !< number of SVF normalization thresholds to report

!-- radiation related arrays to be used in radiation_interaction routine
    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rad_sw_in_dir    !< direct sw radiation
    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rad_sw_in_diff   !< diffusion sw radiation
    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rad_lw_in_diff   !< diffusion lw radiation

!-- parameters required for RRTMG lower boundary condition
    REAL(wp)                   :: albedo_urb      !< albedo value retuned to RRTMG boundary cond.
    REAL(wp)                   :: emissivity_urb  !< emissivity value retuned to RRTMG boundary cond.
    REAL(wp)                   :: t_rad_urb       !< temperature value retuned to RRTMG boundary cond.

!-- type for calculation of svf
    TYPE t_svf
        INTEGER(iwp)                               :: isurflt           !< 
        INTEGER(iwp)                               :: isurfs            !< 
        REAL(wp)                                   :: rsvf              !< 
        REAL(wp)                                   :: rtransp           !< 
    END TYPE

!-- type for calculation of csf
    TYPE t_csf
        INTEGER(iwp)                               :: ip                !< 
        INTEGER(iwp)                               :: itx               !< 
        INTEGER(iwp)                               :: ity               !< 
        INTEGER(iwp)                               :: itz               !< 
        INTEGER(iwp)                               :: isurfs            !< Idx of source face / -1 for sky
        REAL(wp)                                   :: rcvf              !< Canopy view factor for faces /
                                                                        !< canopy sink factor for sky (-1)
    END TYPE

!-- arrays storing the values of USM
    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  svfsurf          !< svfsurf[:,isvf] = index of target and source surface for svf[isvf]
    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  svf              !< array of shape view factors+direct irradiation factors for local surfaces
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins          !< array of sw radiation falling to local surface after i-th reflection
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl          !< array of lw radiation for local surface after i-th reflection

    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  skyvf            !< array of sky view factor for each local surface
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  skyvft           !< array of sky view factor including transparency for each local surface
    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsitrans         !< dsidir[isvfl,i] = path transmittance of i-th
                                                                        !< direction of direct solar irradiance per target surface
    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsitransc        !< dtto per plant canopy box
    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsidir           !< dsidir[:,i] = unit vector of i-th
                                                                        !< direction of direct solar irradiance
    INTEGER(iwp)                                   ::  ndsidir          !< number of apparent solar directions used
    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  dsidir_rev       !< dsidir_rev[ielev,iazim] = i for dsidir or -1 if not present

    INTEGER(iwp)                                   ::  nmrtbl           !< No. of local grid boxes for which MRT is calculated
    INTEGER(iwp)                                   ::  nmrtf            !< number of MRT factors for local processor
    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  mrtbl            !< coordinates of i-th local MRT box - surfl[:,i] = [z, y, x]
    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  mrtfsurf         !< mrtfsurf[:,imrtf] = index of target MRT box and source surface for mrtf[imrtf]
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtf             !< array of MRT factors for each local MRT box
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtft            !< array of MRT factors including transparency for each local MRT box
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtsky           !< array of sky view factor for each local MRT box
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtskyt          !< array of sky view factor including transparency for each local MRT box
    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  mrtdsit          !< array of direct solar transparencies for each local MRT box
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinsw          !< mean SW radiant flux for each MRT box
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinlw          !< mean LW radiant flux for each MRT box
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrt              !< mean radiant temperature for each MRT box
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinsw_av       !< time average mean SW radiant flux for each MRT box
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinlw_av       !< time average mean LW radiant flux for each MRT box
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrt_av           !< time average mean radiant temperature for each MRT box

    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw         !< array of sw radiation falling to local surface including radiation from reflections
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw         !< array of lw radiation falling to local surface including radiation from reflections
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir      !< array of direct sw radiation falling to local surface
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif      !< array of diffuse sw radiation from sky and model boundary falling to local surface
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif      !< array of diffuse lw radiation from sky and model boundary falling to local surface
    
                                                                        !< Outward radiation is only valid for nonvirtual surfaces
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsl        !< array of reflected sw radiation for local surface in i-th reflection
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutll        !< array of reflected + emitted lw radiation for local surface in i-th reflection
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfouts         !< array of reflected sw radiation for all surfaces in i-th reflection
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutl         !< array of reflected + emitted lw radiation for all surfaces in i-th reflection
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlg         !< global array of incoming lw radiation from plant canopy
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw        !< array of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw        !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfemitlwl      !< array of emitted lw radiation for local surface used to calculate effective surface temperature for radiation model

!-- block variables needed for calculation of the plant canopy model inside the urban surface model
    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  csfsurf          !< csfsurf[:,icsf] = index of target surface and csf grid index for csf[icsf]
    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  csf              !< array of plant canopy sink fators + direct irradiation factors (transparency)
    REAL(wp), DIMENSION(:,:,:), POINTER            ::  sub_lad          !< subset of lad_s within urban surface, transformed to plain Z coordinate
    REAL(wp), DIMENSION(:), POINTER                ::  sub_lad_g        !< sub_lad globalized (used to avoid MPI RMA calls in raytracing)
    REAL(wp)                                       ::  prototype_lad    !< prototype leaf area density for computing effective optical depth
    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nzterr, plantt   !< temporary global arrays for raytracing
    INTEGER(iwp)                                   ::  plantt_max

!-- arrays and variables for calculation of svf and csf
    TYPE(t_svf), DIMENSION(:), POINTER             ::  asvf             !< pointer to growing svc array
    TYPE(t_csf), DIMENSION(:), POINTER             ::  acsf             !< pointer to growing csf array
    TYPE(t_svf), DIMENSION(:), POINTER             ::  amrtf            !< pointer to growing mrtf array
    TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  asvf1, asvf2     !< realizations of svf array
    TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET ::  acsf1, acsf2     !< realizations of csf array
    TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  amrtf1, amrtf2   !< realizations of mftf array
    INTEGER(iwp)                                   ::  nsvfla           !< dimmension of array allocated for storage of svf in local processor
    INTEGER(iwp)                                   ::  ncsfla           !< dimmension of array allocated for storage of csf in local processor
    INTEGER(iwp)                                   ::  nmrtfa           !< dimmension of array allocated for storage of mrt
    INTEGER(iwp)                                   ::  msvf, mcsf, mmrtf!< mod for swapping the growing array
    INTEGER(iwp), PARAMETER                        ::  gasize = 100000_iwp  !< initial size of growing arrays
    REAL(wp), PARAMETER                            ::  grow_factor = 1.4_wp !< growth factor of growing arrays
    INTEGER(iwp)                                   ::  nsvfl            !< number of svf for local processor
    INTEGER(iwp)                                   ::  ncsfl            !< no. of csf in local processor
                                                                        !< needed only during calc_svf but must be here because it is
                                                                        !< shared between subroutines calc_svf and raytrace
    INTEGER(iwp), DIMENSION(:,:,:,:), POINTER      ::  gridsurf         !< reverse index of local surfl[d,k,j,i] (for case rad_angular_discretization)
    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE    ::  gridpcbl         !< reverse index of local pcbl[k,j,i]
    INTEGER(iwp), PARAMETER                        ::  nsurf_type_u = 6 !< number of urban surf types (used in gridsurf)

!-- temporary arrays for calculation of csf in raytracing
    INTEGER(iwp)                                   ::  maxboxesg        !< max number of boxes ray can cross in the domain
    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  boxes            !< coordinates of gridboxes being crossed by ray
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  crlens           !< array of crossing lengths of ray for particular grid boxes
    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  lad_ip           !< array of numbers of process where lad is stored 
#if defined( __parallel )
    INTEGER(kind=MPI_ADDRESS_KIND), &
                  DIMENSION(:), ALLOCATABLE        ::  lad_disp         !< array of displaycements of lad in local array of proc lad_ip
    INTEGER(iwp)                                   ::  win_lad          !< MPI RMA window for leaf area density
    INTEGER(iwp)                                   ::  win_gridsurf     !< MPI RMA window for reverse grid surface index
#endif
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  lad_s_ray        !< array of received lad_s for appropriate gridboxes crossed by ray
    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  target_surfl
    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  rt2_track
    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rt2_track_lad
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rt2_track_dist
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rt2_dist

!-- arrays for time averages
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfradnet_av    !< average of net radiation to local surface including radiation from reflections
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw_av      !< average of sw radiation falling to local surface including radiation from reflections
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw_av      !< average of lw radiation falling to local surface including radiation from reflections
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir_av   !< average of direct sw radiation falling to local surface
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif_av   !< average of diffuse sw radiation from sky and model boundary falling to local surface
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif_av   !< average of diffuse lw radiation from sky and model boundary falling to local surface
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswref_av   !< average of sw radiation falling to surface from reflections
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwref_av   !< average of lw radiation falling to surface from reflections
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw_av     !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw_av     !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw_av       !< Average of pcbinlw
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw_av       !< Average of pcbinsw
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir_av    !< Average of pcbinswdir
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif_av    !< Average of pcbinswdif
    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswref_av    !< Average of pcbinswref


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-- Energy balance variables
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-- parameters of the land, roof and wall surfaces
    REAL(wp), DIMENSION(:), ALLOCATABLE            :: albedo_surf        !< albedo of the surface
    REAL(wp), DIMENSION(:), ALLOCATABLE            :: emiss_surf         !< emissivity of the wall surface


    INTERFACE radiation_check_data_output
       MODULE PROCEDURE radiation_check_data_output
    END INTERFACE radiation_check_data_output

    INTERFACE radiation_check_data_output_ts
       MODULE PROCEDURE radiation_check_data_output_ts
    END INTERFACE radiation_check_data_output_ts

    INTERFACE radiation_check_data_output_pr
       MODULE PROCEDURE radiation_check_data_output_pr
    END INTERFACE radiation_check_data_output_pr
  
    INTERFACE radiation_check_parameters
       MODULE PROCEDURE radiation_check_parameters
    END INTERFACE radiation_check_parameters
  
    INTERFACE radiation_clearsky
       MODULE PROCEDURE radiation_clearsky
    END INTERFACE radiation_clearsky
 
    INTERFACE radiation_constant
       MODULE PROCEDURE radiation_constant
    END INTERFACE radiation_constant
 
    INTERFACE radiation_control
       MODULE PROCEDURE radiation_control
    END INTERFACE radiation_control

    INTERFACE radiation_3d_data_averaging
       MODULE PROCEDURE radiation_3d_data_averaging
    END INTERFACE radiation_3d_data_averaging

    INTERFACE radiation_data_output_2d
       MODULE PROCEDURE radiation_data_output_2d
    END INTERFACE radiation_data_output_2d

    INTERFACE radiation_data_output_3d
       MODULE PROCEDURE radiation_data_output_3d
    END INTERFACE radiation_data_output_3d

    INTERFACE radiation_data_output_mask
       MODULE PROCEDURE radiation_data_output_mask
    END INTERFACE radiation_data_output_mask

    INTERFACE radiation_define_netcdf_grid
       MODULE PROCEDURE radiation_define_netcdf_grid
    END INTERFACE radiation_define_netcdf_grid

    INTERFACE radiation_header
       MODULE PROCEDURE radiation_header
    END INTERFACE radiation_header 
  
    INTERFACE radiation_init
       MODULE PROCEDURE radiation_init
    END INTERFACE radiation_init

    INTERFACE radiation_parin
       MODULE PROCEDURE radiation_parin
    END INTERFACE radiation_parin
    
    INTERFACE radiation_rrtmg
       MODULE PROCEDURE radiation_rrtmg
    END INTERFACE radiation_rrtmg

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

    INTERFACE radiation_rrd_local
       MODULE PROCEDURE radiation_rrd_local
    END INTERFACE radiation_rrd_local

    INTERFACE radiation_wrd_local
       MODULE PROCEDURE radiation_wrd_local
    END INTERFACE radiation_wrd_local

    INTERFACE radiation_interaction
       MODULE PROCEDURE radiation_interaction
    END INTERFACE radiation_interaction

    INTERFACE radiation_interaction_init
       MODULE PROCEDURE radiation_interaction_init
    END INTERFACE radiation_interaction_init
 
    INTERFACE radiation_presimulate_solar_pos
       MODULE PROCEDURE radiation_presimulate_solar_pos
    END INTERFACE radiation_presimulate_solar_pos

    INTERFACE radiation_calc_svf
       MODULE PROCEDURE radiation_calc_svf
    END INTERFACE radiation_calc_svf

    INTERFACE radiation_write_svf
       MODULE PROCEDURE radiation_write_svf
    END INTERFACE radiation_write_svf

    INTERFACE radiation_read_svf
       MODULE PROCEDURE radiation_read_svf
    END INTERFACE radiation_read_svf


    SAVE

    PRIVATE

!
!-- Public functions / NEEDS SORTING
    PUBLIC radiation_check_data_output, radiation_check_data_output_pr,        &
           radiation_check_data_output_ts,                                     &
           radiation_check_parameters, radiation_control,                      &
           radiation_header, radiation_init, radiation_parin,                  &
           radiation_3d_data_averaging,                                        &
           radiation_data_output_2d, radiation_data_output_3d,                 &
           radiation_define_netcdf_grid, radiation_wrd_local,                  &
           radiation_rrd_local, radiation_data_output_mask,                    &
           radiation_calc_svf, radiation_write_svf,                            &
           radiation_interaction, radiation_interaction_init,                  &
           radiation_read_svf, radiation_presimulate_solar_pos

    
!
!-- Public variables and constants / NEEDS SORTING
    PUBLIC albedo, albedo_type, decl_1, decl_2, decl_3, dots_rad, dt_radiation,&
           emissivity, force_radiation_call, lat, lon, mrt_geom_human,         &
           mrt_include_sw, mrt_nlevels, mrtbl, mrtinsw, mrtinlw, nmrtbl,       &
           rad_net_av, radiation, radiation_scheme, rad_lw_in,                 &
           rad_lw_in_av, rad_lw_out, rad_lw_out_av,                            &
           rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
           rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, solar_constant,           &
           skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
           cos_zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,   &
           idir, jdir, kdir, id, iz, iy, ix,                                   &
           iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
           iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
           nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf,                 &
           idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
           radiation_interactions, startwall, startland, endland, endwall,     &
           skyvf, skyvft, radiation_interactions_on, average_radiation,        &
           rad_sw_in_diff, rad_sw_in_dir


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

 CONTAINS


!------------------------------------------------------------------------------!
! Description:
! ------------
!> This subroutine controls the calls of the radiation schemes
!------------------------------------------------------------------------------!
    SUBROUTINE radiation_control
 
 
       IMPLICIT NONE


       SELECT CASE ( TRIM( radiation_scheme ) )

          CASE ( 'constant' )
             CALL radiation_constant
          
          CASE ( 'clear-sky' ) 
             CALL radiation_clearsky
       
          CASE ( 'rrtmg' )
             CALL radiation_rrtmg

          CASE DEFAULT

       END SELECT


    END SUBROUTINE radiation_control

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Check data output for radiation model
!------------------------------------------------------------------------------!
    SUBROUTINE radiation_check_data_output( variable, unit, i, ilen, k )
 
 
       USE control_parameters,                                                 &
           ONLY: data_output, message_string

       IMPLICIT NONE

       CHARACTER (LEN=*) ::  unit          !<
       CHARACTER (LEN=*) ::  variable      !<

       INTEGER(iwp) :: i, k
       INTEGER(iwp) :: ilen
       CHARACTER(LEN=varnamelength) :: var          !< TRIM(variable)

       var = TRIM(variable)

!--    first process diractional variables
       IF ( var(1:12) == 'rtm_rad_net_'  .OR.  var(1:13) == 'rtm_rad_insw_'  .OR.        &
            var(1:13) == 'rtm_rad_inlw_'  .OR.  var(1:16) == 'rtm_rad_inswdir_'  .OR.    &
            var(1:16) == 'rtm_rad_inswdif_'  .OR.  var(1:16) == 'rtm_rad_inswref_'  .OR. &
            var(1:16) == 'rtm_rad_inlwdif_'  .OR.  var(1:16) == 'rtm_rad_inlwref_'  .OR. &
            var(1:14) == 'rtm_rad_outsw_'  .OR.  var(1:14) == 'rtm_rad_outlw_'  .OR.     &
            var(1:14) == 'rtm_rad_ressw_'  .OR.  var(1:14) == 'rtm_rad_reslw_'  ) THEN
          IF ( .NOT.  radiation ) THEN
                message_string = 'output of "' // TRIM( var ) // '" require'&
                                 // 's radiation = .TRUE.'
                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
          ENDIF
          unit = 'W/m2'
       ELSE IF ( var(1:7) == 'rtm_svf'  .OR.  var(1:7) == 'rtm_dif'  .OR.                &
                 var(1:9) == 'rtm_skyvf' .OR. var(1:9) == 'rtm_skyvft' )  THEN
          IF ( .NOT.  radiation ) THEN
                message_string = 'output of "' // TRIM( var ) // '" require'&
                                 // 's radiation = .TRUE.'
                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
          ENDIF
          unit = '1'
       ELSE
!--       non-directional variables
          SELECT CASE ( TRIM( var ) )
             CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out', &
                    'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out'  )
                IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
                   message_string = '"output of "' // TRIM( var ) // '" requi' // &
                                    'res radiation = .TRUE. and ' //              &
                                    'radiation_scheme = "rrtmg"'
                   CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
                ENDIF
                unit = 'K/h'

             CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
                    'rrtm_asdir*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',     &
                    'rad_sw_out*')
                IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
                   ! Workaround for masked output (calls with i=ilen=k=0)
                   unit = 'illegal'
                   RETURN
                ENDIF
                IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
                   message_string = 'illegal value for data_output: "' //         &
                                    TRIM( var ) // '" & only 2d-horizontal ' //   &
                                    'cross sections are allowed for this value'
                   CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
                ENDIF
                IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
                   IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
                        TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
                        TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
                        TRIM( var ) == 'rrtm_asdir*'      )                       &
                   THEN
                      message_string = 'output of "' // TRIM( var ) // '" require'&
                                       // 's radiation = .TRUE. and radiation_sch'&
                                       // 'eme = "rrtmg"'
                      CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
                   ENDIF
                ENDIF

                IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'
                IF ( TRIM( var ) == 'rad_lw_in*'    ) unit = 'W/m2'
                IF ( TRIM( var ) == 'rad_lw_out*'   ) unit = 'W/m2'
                IF ( TRIM( var ) == 'rad_sw_in*'    ) unit = 'W/m2'
                IF ( TRIM( var ) == 'rad_sw_out*'   ) unit = 'W/m2'
                IF ( TRIM( var ) == 'rad_sw_in'     ) unit = 'W/m2'
                IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''
                IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = ''
                IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = ''
                IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = ''

             CASE ( 'rtm_rad_pc_inlw', 'rtm_rad_pc_insw', 'rtm_rad_pc_inswdir', &
                    'rtm_rad_pc_inswdif', 'rtm_rad_pc_inswref')
                IF ( .NOT.  radiation ) THEN
                   message_string = 'output of "' // TRIM( var ) // '" require'&
                                    // 's radiation = .TRUE.'
                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
                ENDIF
                unit = 'W'

             CASE ( 'rtm_mrt', 'rtm_mrt_sw', 'rtm_mrt_lw'  )
                IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
                   ! Workaround for masked output (calls with i=ilen=k=0)
                   unit = 'illegal'
                   RETURN
                ENDIF

                IF ( .NOT.  radiation ) THEN
                   message_string = 'output of "' // TRIM( var ) // '" require'&
                                    // 's radiation = .TRUE.'
                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
                ENDIF
                IF ( mrt_nlevels == 0 ) THEN
                   message_string = 'output of "' // TRIM( var ) // '" require'&
                                    // 's mrt_nlevels > 0'
                   CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
                ENDIF
                IF ( TRIM( var ) == 'rtm_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
                   message_string = 'output of "' // TRIM( var ) // '" require'&
                                    // 's rtm_mrt_sw = .TRUE.'
                   CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
                ENDIF
                IF ( TRIM( var ) == 'rtm_mrt' ) THEN
                   unit = 'K'
                ELSE
                   unit = 'W m-2'
                ENDIF

             CASE DEFAULT
                unit = 'illegal'

          END SELECT
       ENDIF

    END SUBROUTINE radiation_check_data_output


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Set module-specific timeseries units and labels
!------------------------------------------------------------------------------!
 SUBROUTINE radiation_check_data_output_ts( dots_max, dots_num )


    INTEGER(iwp),      INTENT(IN)     ::  dots_max
    INTEGER(iwp),      INTENT(INOUT)  ::  dots_num

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

!
!-- Temporary solution to add LSM and radiation time series to the default
!-- output
    IF ( land_surface  .OR.  radiation )  THEN
       IF ( TRIM( radiation_scheme ) == 'rrtmg' )  THEN
          dots_num = dots_num + 15
       ELSE
          dots_num = dots_num + 11
       ENDIF
    ENDIF


 END SUBROUTINE radiation_check_data_output_ts

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Check data output of profiles for radiation model
!------------------------------------------------------------------------------!  
    SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit,      &
               dopr_unit )
 
       USE arrays_3d,                                                          &
           ONLY: zu

       USE control_parameters,                                                 &
           ONLY: data_output_pr, message_string

       USE indices

       USE profil_parameter

       USE statistics

       IMPLICIT NONE
   
       CHARACTER (LEN=*) ::  unit      !< 
       CHARACTER (LEN=*) ::  variable  !< 
       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
 
       INTEGER(iwp) ::  var_count     !< 

       SELECT CASE ( TRIM( variable ) )
       
         CASE ( 'rad_net' )
             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
             THEN
                message_string = 'data_output_pr = ' //                        &
                                 TRIM( data_output_pr(var_count) ) // ' is' // &
                                 'not available for radiation = .FALSE. or ' //&
                                 'radiation_scheme = "constant"'
                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
             ELSE
                dopr_index(var_count) = 99
                dopr_unit  = 'W/m2'
                hom(:,2,99,:)  = SPREAD( zw, 2, statistic_regions+1 )
                unit = dopr_unit
             ENDIF

          CASE ( 'rad_lw_in' )
             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
             THEN
                message_string = 'data_output_pr = ' //                        &
                                 TRIM( data_output_pr(var_count) ) // ' is' // &
                                 'not available for radiation = .FALSE. or ' //&
                                 'radiation_scheme = "constant"'
                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
             ELSE
                dopr_index(var_count) = 100
                dopr_unit  = 'W/m2'
                hom(:,2,100,:)  = SPREAD( zw, 2, statistic_regions+1 )
                unit = dopr_unit 
             ENDIF

          CASE ( 'rad_lw_out' )
             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
             THEN
                message_string = 'data_output_pr = ' //                        &
                                 TRIM( data_output_pr(var_count) ) // ' is' // &
                                 'not available for radiation = .FALSE. or ' //&
                                 'radiation_scheme = "constant"'
                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
             ELSE
                dopr_index(var_count) = 101
                dopr_unit  = 'W/m2'
                hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions+1 )
                unit = dopr_unit   
             ENDIF

          CASE ( 'rad_sw_in' )
             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
             THEN
                message_string = 'data_output_pr = ' //                        &
                                 TRIM( data_output_pr(var_count) ) // ' is' // &
                                 'not available for radiation = .FALSE. or ' //&
                                 'radiation_scheme = "constant"'
                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
             ELSE
                dopr_index(var_count) = 102
                dopr_unit  = 'W/m2'
                hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
                unit = dopr_unit
             ENDIF

          CASE ( 'rad_sw_out')
             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
             THEN
                message_string = 'data_output_pr = ' //                        &
                                 TRIM( data_output_pr(var_count) ) // ' is' // &
                                 'not available for radiation = .FALSE. or ' //&
                                 'radiation_scheme = "constant"'
                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
             ELSE
                dopr_index(var_count) = 103
                dopr_unit  = 'W/m2'
                hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
                unit = dopr_unit
             ENDIF

          CASE ( 'rad_lw_cs_hr' )
             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
             THEN
                message_string = 'data_output_pr = ' //                        &
                                 TRIM( data_output_pr(var_count) ) // ' is' // &
                                 'not available for radiation = .FALSE. or ' //&
                                 'radiation_scheme /= "rrtmg"'
                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
             ELSE
                dopr_index(var_count) = 104
                dopr_unit  = 'K/h'
                hom(:,2,104,:)  = SPREAD( zu, 2, statistic_regions+1 )
                unit = dopr_unit
             ENDIF

          CASE ( 'rad_lw_hr' )
             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
             THEN
                message_string = 'data_output_pr = ' //                        &
                                 TRIM( data_output_pr(var_count) ) // ' is' // &
                                 'not available for radiation = .FALSE. or ' //&
                                 'radiation_scheme /= "rrtmg"'
                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
             ELSE
                dopr_index(var_count) = 105
                dopr_unit  = 'K/h'
                hom(:,2,105,:)  = SPREAD( zu, 2, statistic_regions+1 )
                unit = dopr_unit
             ENDIF

          CASE ( 'rad_sw_cs_hr' )
             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
             THEN
                message_string = 'data_output_pr = ' //                        &
                                 TRIM( data_output_pr(var_count) ) // ' is' // &
                                 'not available for radiation = .FALSE. or ' //&
                                 'radiation_scheme /= "rrtmg"'
                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
             ELSE
                dopr_index(var_count) = 106
                dopr_unit  = 'K/h'
                hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions+1 )
                unit = dopr_unit
             ENDIF

          CASE ( 'rad_sw_hr' )
             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
             THEN
                message_string = 'data_output_pr = ' //                        &
                                 TRIM( data_output_pr(var_count) ) // ' is' // &
                                 'not available for radiation = .FALSE. or ' //&
                                 'radiation_scheme /= "rrtmg"'
                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
             ELSE
                dopr_index(var_count) = 107
                dopr_unit  = 'K/h'
                hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions+1 )
                unit = dopr_unit
             ENDIF


          CASE DEFAULT
             unit = 'illegal'

       END SELECT


    END SUBROUTINE radiation_check_data_output_pr
 
 
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Check parameters routine for radiation model
!------------------------------------------------------------------------------!
    SUBROUTINE radiation_check_parameters

       USE control_parameters,                                                 &
           ONLY: land_surface, message_string, urban_surface

       USE netcdf_data_input_mod,                                              &
           ONLY:  input_pids_static                 
    
       IMPLICIT NONE
       
!
!--    In case no urban-surface or land-surface model is applied, usage of 
!--    a radiation model make no sense.         
       IF ( .NOT. land_surface  .AND.  .NOT. urban_surface )  THEN
          message_string = 'Usage of radiation module is only allowed if ' //  &
                           'land-surface and/or urban-surface model is applied.'
          CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
       ENDIF

       IF ( radiation_scheme /= 'constant'   .AND.                             &
            radiation_scheme /= 'clear-sky'  .AND.                             &
            radiation_scheme /= 'rrtmg' )  THEN
          message_string = 'unknown radiation_scheme = '//                     &
                           TRIM( radiation_scheme )
          CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
#if ! defined ( __rrtmg )
          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
                           'compilation of PALM with pre-processor ' //        &
                           'directive -D__rrtmg'
          CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
#endif
#if defined ( __rrtmg ) && ! defined( __netcdf )
          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
                           'the use of NetCDF (preprocessor directive ' //     &
                           '-D__netcdf'
          CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
#endif

       ENDIF
!
!--    Checks performed only if data is given via namelist only. 
       IF ( .NOT. input_pids_static )  THEN
          IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.          &
               radiation_scheme == 'clear-sky')  THEN
             message_string = 'radiation_scheme = "clear-sky" in combination'//& 
                              'with albedo_type = 0 requires setting of'//     &
                              'albedo /= 9999999.9'
             CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
          ENDIF

          IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.     &
             ( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
          .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& 
             ) ) THEN
             message_string = 'radiation_scheme = "rrtmg" in combination' //   & 
                              'with albedo_type = 0 requires setting of ' //   &
                              'albedo_lw_dif /= 9999999.9' //                  &
                              'albedo_lw_dir /= 9999999.9' //                  &
                              'albedo_sw_dif /= 9999999.9 and' //              &
                              'albedo_sw_dir /= 9999999.9'
             CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
          ENDIF
       ENDIF
!
!--    Parallel rad_angular_discretization without raytrace_mpi_rma is not implemented
#if defined( __parallel )      
       IF ( rad_angular_discretization  .AND.  .NOT. raytrace_mpi_rma )  THEN
          message_string = 'rad_angular_discretization can only be used ' //  &
                           'together with raytrace_mpi_rma or when ' //  &
                           'no parallelization is applied.'
          CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
       ENDIF
#endif

       IF ( cloud_droplets  .AND.   radiation_scheme == 'rrtmg'  .AND.         &
            average_radiation ) THEN
          message_string = 'average_radiation = .T. with radiation_scheme'//   & 
                           '= "rrtmg" in combination cloud_droplets = .T.'//   &
                           'is not implementd'
          CALL message( 'check_parameters', 'PA0560', 1, 2, 0, 6, 0 )
       ENDIF

!
!--    Incialize svf normalization reporting histogram
       svfnorm_report_num = 1
       DO WHILE ( svfnorm_report_thresh(svfnorm_report_num) < 1e20_wp          &
                   .AND. svfnorm_report_num <= 30 )
          svfnorm_report_num = svfnorm_report_num + 1
       ENDDO
       svfnorm_report_num = svfnorm_report_num - 1
!
!--    Check for dt_radiation
       IF ( dt_radiation <= 0.0 )  THEN
          message_string = 'dt_radiation must be > 0.0'  
          CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 ) 
       ENDIF
 
    END SUBROUTINE radiation_check_parameters 
 
 
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Initialization of the radiation model
!------------------------------------------------------------------------------!
    SUBROUTINE radiation_init
    
       IMPLICIT NONE

       INTEGER(iwp) ::  i         !< running index x-direction 
       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) ::  l         !< running index for orientation of vertical surfaces 
       INTEGER(iwp) ::  m         !< running index for surface elements
#if defined( __rrtmg )
       INTEGER(iwp) ::  ind_type  !< running index for subgrid-surface tiles
#endif

!
!--    Activate radiation_interactions according to the existence of vertical surfaces and/or trees.
!--    The namelist parameter radiation_interactions_on can override this behavior.
!--    (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in 
!--    init_surface_arrays.)
       IF ( radiation_interactions_on )  THEN
          IF ( vertical_surfaces_exist  .OR.  plant_canopy )  THEN
             radiation_interactions    = .TRUE.
             average_radiation         = .TRUE.
          ELSE
             radiation_interactions_on = .FALSE.   !< reset namelist parameter: no interactions 
                                                   !< calculations necessary in case of flat surface
          ENDIF
       ELSEIF ( vertical_surfaces_exist  .OR.  plant_canopy )  THEN
          message_string = 'radiation_interactions_on is set to .FALSE. although '     // &
                           'vertical surfaces and/or trees exist. The model will run ' // &
                           'without RTM (no shadows, no radiation reflections)'
          CALL message( 'init_3d_model', 'PA0348', 0, 1, 0, 6, 0 )
       ENDIF
!
!--    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

!
!--    Initialize radiation model
       CALL location_message( 'initializing radiation model', .FALSE. )

!
!--    Allocate array for storing the surface net radiation
       IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_net )  .AND.                      &
                  surf_lsm_h%ns > 0  )   THEN
          ALLOCATE( surf_lsm_h%rad_net(1:surf_lsm_h%ns) )
          surf_lsm_h%rad_net = 0.0_wp 
       ENDIF
       IF ( .NOT. ALLOCATED ( surf_usm_h%rad_net )  .AND.                      &
                  surf_usm_h%ns > 0  )  THEN
          ALLOCATE( surf_usm_h%rad_net(1:surf_usm_h%ns) )
          surf_usm_h%rad_net = 0.0_wp 
       ENDIF
       DO  l = 0, 3
          IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_net )  .AND.                &
                     surf_lsm_v(l)%ns > 0  )  THEN
             ALLOCATE( surf_lsm_v(l)%rad_net(1:surf_lsm_v(l)%ns) )
             surf_lsm_v(l)%rad_net = 0.0_wp 
          ENDIF
          IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_net )  .AND.                &
                     surf_usm_v(l)%ns > 0  )  THEN
             ALLOCATE( surf_usm_v(l)%rad_net(1:surf_usm_v(l)%ns) )
             surf_usm_v(l)%rad_net = 0.0_wp 
          ENDIF
       ENDDO


!
!--    Allocate array for storing the surface longwave (out) radiation change
       IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_lw_out_change_0 )  .AND.          &
                  surf_lsm_h%ns > 0  )   THEN
          ALLOCATE( surf_lsm_h%rad_lw_out_change_0(1:surf_lsm_h%ns) )
          surf_lsm_h%rad_lw_out_change_0 = 0.0_wp 
       ENDIF
       IF ( .NOT. ALLOCATED ( surf_usm_h%rad_lw_out_change_0 )  .AND.          &
                  surf_usm_h%ns > 0  )  THEN
          ALLOCATE( surf_usm_h%rad_lw_out_change_0(1:surf_usm_h%ns) )
          surf_usm_h%rad_lw_out_change_0 = 0.0_wp 
       ENDIF
       DO  l = 0, 3
          IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_lw_out_change_0 )  .AND.    &
                     surf_lsm_v(l)%ns > 0  )  THEN
             ALLOCATE( surf_lsm_v(l)%rad_lw_out_change_0(1:surf_lsm_v(l)%ns) )
             surf_lsm_v(l)%rad_lw_out_change_0 = 0.0_wp 
          ENDIF
          IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_lw_out_change_0 )  .AND.    &
                     surf_usm_v(l)%ns > 0  )  THEN
             ALLOCATE( surf_usm_v(l)%rad_lw_out_change_0(1:surf_usm_v(l)%ns) )
             surf_usm_v(l)%rad_lw_out_change_0 = 0.0_wp 
          ENDIF
       ENDDO

!
!--    Allocate surface arrays for incoming/outgoing short/longwave radiation
       IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_sw_in )  .AND.                    &
                  surf_lsm_h%ns > 0  )   THEN
          ALLOCATE( surf_lsm_h%rad_sw_in(1:surf_lsm_h%ns)  )
          ALLOCATE( surf_lsm_h%rad_sw_out(1:surf_lsm_h%ns) )
          ALLOCATE( surf_lsm_h%rad_sw_dir(1:surf_lsm_h%ns) )
          ALLOCATE( surf_lsm_h%rad_sw_dif(1:surf_lsm_h%ns) )
          ALLOCATE( surf_lsm_h%rad_sw_ref(1:surf_lsm_h%ns) )
          ALLOCATE( surf_lsm_h%rad_sw_res(1:surf_lsm_h%ns) )
          ALLOCATE( surf_lsm_h%rad_lw_in(1:surf_lsm_h%ns)  )
          ALLOCATE( surf_lsm_h%rad_lw_out(1:surf_lsm_h%ns) )
          ALLOCATE( surf_lsm_h%rad_lw_dif(1:surf_lsm_h%ns) )
          ALLOCATE( surf_lsm_h%rad_lw_ref(1:surf_lsm_h%ns) )
          ALLOCATE( surf_lsm_h%rad_lw_res(1:surf_lsm_h%ns) )
          surf_lsm_h%rad_sw_in  = 0.0_wp 
          surf_lsm_h%rad_sw_out = 0.0_wp 
          surf_lsm_h%rad_sw_dir = 0.0_wp 
          surf_lsm_h%rad_sw_dif = 0.0_wp 
          surf_lsm_h%rad_sw_ref = 0.0_wp 
          surf_lsm_h%rad_sw_res = 0.0_wp 
          surf_lsm_h%rad_lw_in  = 0.0_wp 
          surf_lsm_h%rad_lw_out = 0.0_wp 
          surf_lsm_h%rad_lw_dif = 0.0_wp 
          surf_lsm_h%rad_lw_ref = 0.0_wp 
          surf_lsm_h%rad_lw_res = 0.0_wp 
       ENDIF
       IF ( .NOT. ALLOCATED ( surf_usm_h%rad_sw_in )  .AND.                    &
                  surf_usm_h%ns > 0  )  THEN
          ALLOCATE( surf_usm_h%rad_sw_in(1:surf_usm_h%ns)  )
          ALLOCATE( surf_usm_h%rad_sw_out(1:surf_usm_h%ns) )
          ALLOCATE( surf_usm_h%rad_sw_dir(1:surf_usm_h%ns) )
          ALLOCATE( surf_usm_h%rad_sw_dif(1:surf_usm_h%ns) )
          ALLOCATE( surf_usm_h%rad_sw_ref(1:surf_usm_h%ns) )
          ALLOCATE( surf_usm_h%rad_sw_res(1:surf_usm_h%ns) )
          ALLOCATE( surf_usm_h%rad_lw_in(1:surf_usm_h%ns)  )
          ALLOCATE( surf_usm_h%rad_lw_out(1:surf_usm_h%ns) )
          ALLOCATE( surf_usm_h%rad_lw_dif(1:surf_usm_h%ns) )
          ALLOCATE( surf_usm_h%rad_lw_ref(1:surf_usm_h%ns) )
          ALLOCATE( surf_usm_h%rad_lw_res(1:surf_usm_h%ns) )
          surf_usm_h%rad_sw_in  = 0.0_wp 
          surf_usm_h%rad_sw_out = 0.0_wp 
          surf_usm_h%rad_sw_dir = 0.0_wp 
          surf_usm_h%rad_sw_dif = 0.0_wp 
          surf_usm_h%rad_sw_ref = 0.0_wp 
          surf_usm_h%rad_sw_res = 0.0_wp 
          surf_usm_h%rad_lw_in  = 0.0_wp 
          surf_usm_h%rad_lw_out = 0.0_wp 
          surf_usm_h%rad_lw_dif = 0.0_wp 
          surf_usm_h%rad_lw_ref = 0.0_wp 
          surf_usm_h%rad_lw_res = 0.0_wp 
       ENDIF
       DO  l = 0, 3
          IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_sw_in )  .AND.              &
                     surf_lsm_v(l)%ns > 0  )  THEN
             ALLOCATE( surf_lsm_v(l)%rad_sw_in(1:surf_lsm_v(l)%ns)  )
             ALLOCATE( surf_lsm_v(l)%rad_sw_out(1:surf_lsm_v(l)%ns) )
             ALLOCATE( surf_lsm_v(l)%rad_sw_dir(1:surf_lsm_v(l)%ns) )
             ALLOCATE( surf_lsm_v(l)%rad_sw_dif(1:surf_lsm_v(l)%ns) )
             ALLOCATE( surf_lsm_v(l)%rad_sw_ref(1:surf_lsm_v(l)%ns) )
             ALLOCATE( surf_lsm_v(l)%rad_sw_res(1:surf_lsm_v(l)%ns) )

             ALLOCATE( surf_lsm_v(l)%rad_lw_in(1:surf_lsm_v(l)%ns)  )
             ALLOCATE( surf_lsm_v(l)%rad_lw_out(1:surf_lsm_v(l)%ns) )
             ALLOCATE( surf_lsm_v(l)%rad_lw_dif(1:surf_lsm_v(l)%ns) )
             ALLOCATE( surf_lsm_v(l)%rad_lw_ref(1:surf_lsm_v(l)%ns) )
             ALLOCATE( surf_lsm_v(l)%rad_lw_res(1:surf_lsm_v(l)%ns) )

             surf_lsm_v(l)%rad_sw_in  = 0.0_wp 
             surf_lsm_v(l)%rad_sw_out = 0.0_wp
             surf_lsm_v(l)%rad_sw_dir = 0.0_wp
             surf_lsm_v(l)%rad_sw_dif = 0.0_wp
             surf_lsm_v(l)%rad_sw_ref = 0.0_wp
             surf_lsm_v(l)%rad_sw_res = 0.0_wp

             surf_lsm_v(l)%rad_lw_in  = 0.0_wp 
             surf_lsm_v(l)%rad_lw_out = 0.0_wp 
             surf_lsm_v(l)%rad_lw_dif = 0.0_wp 
             surf_lsm_v(l)%rad_lw_ref = 0.0_wp 
             surf_lsm_v(l)%rad_lw_res = 0.0_wp 
          ENDIF
          IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_sw_in )  .AND.              &
                     surf_usm_v(l)%ns > 0  )  THEN
             ALLOCATE( surf_usm_v(l)%rad_sw_in(1:surf_usm_v(l)%ns)  )
             ALLOCATE( surf_usm_v(l)%rad_sw_out(1:surf_usm_v(l)%ns) )
             ALLOCATE( surf_usm_v(l)%rad_sw_dir(1:surf_usm_v(l)%ns) )
             ALLOCATE( surf_usm_v(l)%rad_sw_dif(1:surf_usm_v(l)%ns) )
             ALLOCATE( surf_usm_v(l)%rad_sw_ref(1:surf_usm_v(l)%ns) )
             ALLOCATE( surf_usm_v(l)%rad_sw_res(1:surf_usm_v(l)%ns) )
             ALLOCATE( surf_usm_v(l)%rad_lw_in(1:surf_usm_v(l)%ns)  )
             ALLOCATE( surf_usm_v(l)%rad_lw_out(1:surf_usm_v(l)%ns) )
             ALLOCATE( surf_usm_v(l)%rad_lw_dif(1:surf_usm_v(l)%ns) )
             ALLOCATE( surf_usm_v(l)%rad_lw_ref(1:surf_usm_v(l)%ns) )
             ALLOCATE( surf_usm_v(l)%rad_lw_res(1:surf_usm_v(l)%ns) )
             surf_usm_v(l)%rad_sw_in  = 0.0_wp 
             surf_usm_v(l)%rad_sw_out = 0.0_wp
             surf_usm_v(l)%rad_sw_dir = 0.0_wp
             surf_usm_v(l)%rad_sw_dif = 0.0_wp
             surf_usm_v(l)%rad_sw_ref = 0.0_wp
             surf_usm_v(l)%rad_sw_res = 0.0_wp
             surf_usm_v(l)%rad_lw_in  = 0.0_wp 
             surf_usm_v(l)%rad_lw_out = 0.0_wp 
             surf_usm_v(l)%rad_lw_dif = 0.0_wp 
             surf_usm_v(l)%rad_lw_ref = 0.0_wp 
             surf_usm_v(l)%rad_lw_res = 0.0_wp 
          ENDIF
       ENDDO
!
!--    Fix net radiation in case of radiation_scheme = 'constant'
       IF ( radiation_scheme == 'constant' )  THEN
          IF ( ALLOCATED( surf_lsm_h%rad_net ) )                               &
             surf_lsm_h%rad_net    = net_radiation
          IF ( ALLOCATED( surf_usm_h%rad_net ) )                               &
             surf_usm_h%rad_net    = net_radiation
!
!--       Todo: weight with inclination angle
          DO  l = 0, 3
             IF ( ALLOCATED( surf_lsm_v(l)%rad_net ) )                         &
                surf_lsm_v(l)%rad_net = net_radiation
             IF ( ALLOCATED( surf_usm_v(l)%rad_net ) )                         &
                surf_usm_v(l)%rad_net = net_radiation
          ENDDO
!          radiation = .FALSE.
!
!--    Calculate orbital constants
       ELSE
          decl_1 = SIN(23.45_wp * pi / 180.0_wp)
          decl_2 = 2.0_wp * pi / 365.0_wp
          decl_3 = decl_2 * 81.0_wp
          lat    = latitude * pi / 180.0_wp
          lon    = longitude * pi / 180.0_wp
       ENDIF

       IF ( radiation_scheme == 'clear-sky'  .OR.                              &
            radiation_scheme == 'constant')  THEN


!
!--       Allocate arrays for incoming/outgoing short/longwave radiation
          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
             ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
          ENDIF
          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
             ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
             ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
          ENDIF
          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
          ENDIF

!
!--       Allocate average arrays for incoming/outgoing short/longwave radiation
          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
             ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
          ENDIF
          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
             ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
          ENDIF

          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
             ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
          ENDIF
          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
             ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
          ENDIF
!
!--       Allocate arrays for broadband albedo, and level 1 initialization
!--       via namelist paramter, unless not already allocated.
          IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )  THEN
             ALLOCATE( surf_lsm_h%albedo(0:2,1:surf_lsm_h%ns)     )
             surf_lsm_h%albedo    = albedo
          ENDIF
          IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )  THEN
             ALLOCATE( surf_usm_h%albedo(0:2,1:surf_usm_h%ns)     )
             surf_usm_h%albedo    = albedo
          ENDIF

          DO  l = 0, 3
             IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )  THEN
                ALLOCATE( surf_lsm_v(l)%albedo(0:2,1:surf_lsm_v(l)%ns) )
                surf_lsm_v(l)%albedo = albedo
             ENDIF
             IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )  THEN
                ALLOCATE( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns) )
                surf_usm_v(l)%albedo = albedo
             ENDIF
          ENDDO
!
!--       Level 2 initialization of broadband albedo via given albedo_type. 
!--       Only if albedo_type is non-zero. In case of urban surface and 
!--       input data is read from ASCII file, albedo_type will be zero, so that
!--       albedo won't be overwritten. 
          DO  m = 1, surf_lsm_h%ns
             IF ( surf_lsm_h%albedo_type(ind_veg_wall,m) /= 0 )                &
                surf_lsm_h%albedo(ind_veg_wall,m) =                            &
                           albedo_pars(2,surf_lsm_h%albedo_type(ind_veg_wall,m))
             IF ( surf_lsm_h%albedo_type(ind_pav_green,m) /= 0 )               &
                surf_lsm_h%albedo(ind_pav_green,m) =                           &
                           albedo_pars(2,surf_lsm_h%albedo_type(ind_pav_green,m))
             IF ( surf_lsm_h%albedo_type(ind_wat_win,m) /= 0 )                 &
                surf_lsm_h%albedo(ind_wat_win,m) =                             &
                           albedo_pars(2,surf_lsm_h%albedo_type(ind_wat_win,m))
          ENDDO
          DO  m = 1, surf_usm_h%ns
             IF ( surf_usm_h%albedo_type(ind_veg_wall,m) /= 0 )                &
                surf_usm_h%albedo(ind_veg_wall,m) =                            &
                           albedo_pars(2,surf_usm_h%albedo_type(ind_veg_wall,m))
             IF ( surf_usm_h%albedo_type(ind_pav_green,m) /= 0 )               &
                surf_usm_h%albedo(ind_pav_green,m) =                           &
                           albedo_pars(2,surf_usm_h%albedo_type(ind_pav_green,m))
             IF ( surf_usm_h%albedo_type(ind_wat_win,m) /= 0 )                 &
                surf_usm_h%albedo(ind_wat_win,m) =                             &
                           albedo_pars(2,surf_usm_h%albedo_type(ind_wat_win,m))
          ENDDO

          DO  l = 0, 3
             DO  m = 1, surf_lsm_v(l)%ns
                IF ( surf_lsm_v(l)%albedo_type(ind_veg_wall,m) /= 0 )          &
                   surf_lsm_v(l)%albedo(ind_veg_wall,m) =                      &
                        albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_veg_wall,m))
                IF ( surf_lsm_v(l)%albedo_type(ind_pav_green,m) /= 0 )         &
                   surf_lsm_v(l)%albedo(ind_pav_green,m) =                     &
                        albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_pav_green,m))
                IF ( surf_lsm_v(l)%albedo_type(ind_wat_win,m) /= 0 )           &
                   surf_lsm_v(l)%albedo(ind_wat_win,m) =                       &
                        albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_wat_win,m))
             ENDDO
             DO  m = 1, surf_usm_v(l)%ns
                IF ( surf_usm_v(l)%albedo_type(ind_veg_wall,m) /= 0 )          &
                   surf_usm_v(l)%albedo(ind_veg_wall,m) =                      &
                        albedo_pars(2,surf_usm_v(l)%albedo_type(ind_veg_wall,m))
                IF ( surf_usm_v(l)%albedo_type(ind_pav_green,m) /= 0 )         &
                   surf_usm_v(l)%albedo(ind_pav_green,m) =                     &
                        albedo_pars(2,surf_usm_v(l)%albedo_type(ind_pav_green,m))
                IF ( surf_usm_v(l)%albedo_type(ind_wat_win,m) /= 0 )           &
                   surf_usm_v(l)%albedo(ind_wat_win,m) =                       &
                        albedo_pars(2,surf_usm_v(l)%albedo_type(ind_wat_win,m))
             ENDDO
          ENDDO

!
!--       Level 3 initialization at grid points where albedo type is zero.
!--       This case, albedo is taken from file. In case of constant radiation
!--       or clear sky, only broadband albedo is given.
          IF ( albedo_pars_f%from_file )  THEN
!
!--          Horizontal surfaces
             DO  m = 1, surf_lsm_h%ns
                i = surf_lsm_h%i(m)
                j = surf_lsm_h%j(m)
                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
                   IF ( surf_lsm_h%albedo_type(ind_veg_wall,m) == 0 )          &
                      surf_lsm_h%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
                   IF ( surf_lsm_h%albedo_type(ind_pav_green,m) == 0 )         &
                      surf_lsm_h%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
                   IF ( surf_lsm_h%albedo_type(ind_wat_win,m) == 0 )           &
                      surf_lsm_h%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
                ENDIF
             ENDDO
             DO  m = 1, surf_usm_h%ns
                i = surf_usm_h%i(m)
                j = surf_usm_h%j(m)
                IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
                   IF ( surf_usm_h%albedo_type(ind_veg_wall,m) == 0 )          &
                      surf_usm_h%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
                   IF ( surf_usm_h%albedo_type(ind_pav_green,m) == 0 )         &
                      surf_usm_h%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
                   IF ( surf_usm_h%albedo_type(ind_wat_win,m) == 0 )           &
                      surf_usm_h%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
                ENDIF
             ENDDO 
!
!--          Vertical surfaces            
             DO  l = 0, 3

                ioff = surf_lsm_v(l)%ioff
                joff = surf_lsm_v(l)%joff
                DO  m = 1, surf_lsm_v(l)%ns
                   i = surf_lsm_v(l)%i(m) + ioff
                   j = surf_lsm_v(l)%j(m) + joff
                   IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
                      IF ( surf_lsm_v(l)%albedo_type(ind_veg_wall,m) == 0 )    &
                         surf_lsm_v(l)%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
                      IF ( surf_lsm_v(l)%albedo_type(ind_pav_green,m) == 0 )   &
                         surf_lsm_v(l)%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
                      IF ( surf_lsm_v(l)%albedo_type(ind_wat_win,m) == 0 )     &
                         surf_lsm_v(l)%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
                   ENDIF
                ENDDO

                ioff = surf_usm_v(l)%ioff
                joff = surf_usm_v(l)%joff
                DO  m = 1, surf_usm_h%ns
                   i = surf_usm_h%i(m) + joff
                   j = surf_usm_h%j(m) + joff
                   IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill )  THEN
                      IF ( surf_usm_v(l)%albedo_type(ind_veg_wall,m) == 0 )    &
                         surf_usm_v(l)%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
                      IF ( surf_usm_v(l)%albedo_type(ind_pav_green,m) == 0 )   &
                         surf_usm_v(l)%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
                      IF ( surf_usm_v(l)%albedo_type(ind_wat_win,m) == 0 )     &
                         surf_lsm_v(l)%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
                   ENDIF
                ENDDO
             ENDDO

          ENDIF  
!
!--    Initialization actions for RRTMG
       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
#if defined ( __rrtmg )
!
!--       Allocate albedos for short/longwave radiation, horizontal surfaces
!--       for wall/green/window (USM) or vegetation/pavement/water surfaces
!--       (LSM). 
          ALLOCATE ( surf_lsm_h%aldif(0:2,1:surf_lsm_h%ns)       )
          ALLOCATE ( surf_lsm_h%aldir(0:2,1:surf_lsm_h%ns)       )
          ALLOCATE ( surf_lsm_h%asdif(0:2,1:surf_lsm_h%ns)       )
          ALLOCATE ( surf_lsm_h%asdir(0:2,1:surf_lsm_h%ns)       )
          ALLOCATE ( surf_lsm_h%rrtm_aldif(0:2,1:surf_lsm_h%ns)  )
          ALLOCATE ( surf_lsm_h%rrtm_aldir(0:2,1:surf_lsm_h%ns)  )
          ALLOCATE ( surf_lsm_h%rrtm_asdif(0:2,1:surf_lsm_h%ns)  )
          ALLOCATE ( surf_lsm_h%rrtm_asdir(0:2,1:surf_lsm_h%ns)  )

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

!
!--       Allocate broadband albedo (temporary for the current radiation 
!--       implementations)
          IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )                            &
             ALLOCATE( surf_lsm_h%albedo(0:2,1:surf_lsm_h%ns)     )
          IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )                            &
             ALLOCATE( surf_usm_h%albedo(0:2,1:surf_usm_h%ns)     )

!
!--       Allocate albedos for short/longwave radiation, vertical surfaces 
          DO  l = 0, 3

             ALLOCATE ( surf_lsm_v(l)%aldif(0:2,1:surf_lsm_v(l)%ns)      )
             ALLOCATE ( surf_lsm_v(l)%aldir(0:2,1:surf_lsm_v(l)%ns)      )
             ALLOCATE ( surf_lsm_v(l)%asdif(0:2,1:surf_lsm_v(l)%ns)      )
             ALLOCATE ( surf_lsm_v(l)%asdir(0:2,1:surf_lsm_v(l)%ns)      )

             ALLOCATE ( surf_lsm_v(l)%rrtm_aldif(0:2,1:surf_lsm_v(l)%ns) )
             ALLOCATE ( surf_lsm_v(l)%rrtm_aldir(0:2,1:surf_lsm_v(l)%ns) )
             ALLOCATE ( surf_lsm_v(l)%rrtm_asdif(0:2,1:surf_lsm_v(l)%ns) )
             ALLOCATE ( surf_lsm_v(l)%rrtm_asdir(0:2,1:surf_lsm_v(l)%ns) )

             ALLOCATE ( surf_usm_v(l)%aldif(0:2,1:surf_usm_v(l)%ns)      )
             ALLOCATE ( surf_usm_v(l)%aldir(0:2,1:surf_usm_v(l)%ns)      )
             ALLOCATE ( surf_usm_v(l)%asdif(0:2,1:surf_usm_v(l)%ns)      )
             ALLOCATE ( surf_usm_v(l)%asdir(0:2,1:surf_usm_v(l)%ns)      )

             ALLOCATE ( surf_usm_v(l)%rrtm_aldif(0:2,1:surf_usm_v(l)%ns) )
             ALLOCATE ( surf_usm_v(l)%rrtm_aldir(0:2,1:surf_usm_v(l)%ns) )
             ALLOCATE ( surf_usm_v(l)%rrtm_asdif(0:2,1:surf_usm_v(l)%ns) )
             ALLOCATE ( surf_usm_v(l)%rrtm_asdir(0:2,1:surf_usm_v(l)%ns) )
!
!--          Allocate broadband albedo (temporary for the current radiation
!--          implementations)
             IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )                    &
                ALLOCATE( surf_lsm_v(l)%albedo(0:2,1:surf_lsm_v(l)%ns) )
             IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )                    &
                ALLOCATE( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns) )

          ENDDO
!
!--       Level 1 initialization of spectral albedos via namelist 
!--       paramters. Please note, this case all surface tiles are initialized 
!--       the same. 
          IF ( surf_lsm_h%ns > 0 )  THEN
             surf_lsm_h%aldif  = albedo_lw_dif
             surf_lsm_h%aldir  = albedo_lw_dir
             surf_lsm_h%asdif  = albedo_sw_dif
             surf_lsm_h%asdir  = albedo_sw_dir
             surf_lsm_h%albedo = albedo_sw_dif
          ENDIF
          IF ( surf_usm_h%ns > 0 )  THEN
             IF ( surf_usm_h%albedo_from_ascii )  THEN
                surf_usm_h%aldif  = surf_usm_h%albedo
                surf_usm_h%aldir  = surf_usm_h%albedo
                surf_usm_h%asdif  = surf_usm_h%albedo
                surf_usm_h%asdir  = surf_usm_h%albedo
             ELSE
                surf_usm_h%aldif  = albedo_lw_dif
                surf_usm_h%aldir  = albedo_lw_dir
                surf_usm_h%asdif  = albedo_sw_dif
                surf_usm_h%asdir  = albedo_sw_dir
                surf_usm_h%albedo = albedo_sw_dif
             ENDIF
          ENDIF

          DO  l = 0, 3

             IF ( surf_lsm_v(l)%ns > 0 )  THEN
                surf_lsm_v(l)%aldif  = albedo_lw_dif
                surf_lsm_v(l)%aldir  = albedo_lw_dir
                surf_lsm_v(l)%asdif  = albedo_sw_dif
                surf_lsm_v(l)%asdir  = albedo_sw_dir
                surf_lsm_v(l)%albedo = albedo_sw_dif
             ENDIF

             IF ( surf_usm_v(l)%ns > 0 )  THEN
                IF ( surf_usm_v(l)%albedo_from_ascii )  THEN
                   surf_usm_v(l)%aldif  = surf_usm_v(l)%albedo
                   surf_usm_v(l)%aldir  = surf_usm_v(l)%albedo
                   surf_usm_v(l)%asdif  = surf_usm_v(l)%albedo
                   surf_usm_v(l)%asdir  = surf_usm_v(l)%albedo
                ELSE
                   surf_usm_v(l)%aldif  = albedo_lw_dif
                   surf_usm_v(l)%aldir  = albedo_lw_dir
                   surf_usm_v(l)%asdif  = albedo_sw_dif
                   surf_usm_v(l)%asdir  = albedo_sw_dir
                ENDIF
             ENDIF
          ENDDO

!
!--       Level 2 initialization of spectral albedos via albedo_type. 
!--       Please note, for natural- and urban-type surfaces, a tile approach 
!--       is applied so that the resulting albedo is calculated via the weighted 
!--       average of respective surface fractions. 
          DO  m = 1, surf_lsm_h%ns
!
!--          Spectral albedos for vegetation/pavement/water surfaces
             DO  ind_type = 0, 2
                IF ( surf_lsm_h%albedo_type(ind_type,m) /= 0 )  THEN
                   surf_lsm_h%aldif(ind_type,m) =                              &
                               albedo_pars(0,surf_lsm_h%albedo_type(ind_type,m))
                   surf_lsm_h%asdif(ind_type,m) =                              &
                               albedo_pars(1,surf_lsm_h%albedo_type(ind_type,m))
                   surf_lsm_h%aldir(ind_type,m) =                              &
                               albedo_pars(0,surf_lsm_h%albedo_type(ind_type,m))
                   surf_lsm_h%asdir(ind_type,m) =                              &
                               albedo_pars(1,surf_lsm_h%albedo_type(ind_type,m))
                   surf_lsm_h%albedo(ind_type,m) =                             &
                               albedo_pars(2,surf_lsm_h%albedo_type(ind_type,m))
                ENDIF
             ENDDO

          ENDDO
!
!--       For urban surface only if albedo has not been already initialized 
!--       in the urban-surface model via the ASCII file.
          IF ( .NOT. surf_usm_h%albedo_from_ascii )  THEN
             DO  m = 1, surf_usm_h%ns
!
!--             Spectral albedos for wall/green/window surfaces
                DO  ind_type = 0, 2
                   IF ( surf_usm_h%albedo_type(ind_type,m) /= 0 )  THEN
                      surf_usm_h%aldif(ind_type,m) =                           &
                               albedo_pars(0,surf_usm_h%albedo_type(ind_type,m))
                      surf_usm_h%asdif(ind_type,m) =                           &
                               albedo_pars(1,surf_usm_h%albedo_type(ind_type,m))
                      surf_usm_h%aldir(ind_type,m) =                           &
                               albedo_pars(0,surf_usm_h%albedo_type(ind_type,m))
                      surf_usm_h%asdir(ind_type,m) =                           &
                               albedo_pars(1,surf_usm_h%albedo_type(ind_type,m))
                      surf_usm_h%albedo(ind_type,m) =                          &
                               albedo_pars(2,surf_usm_h%albedo_type(ind_type,m))
                   ENDIF
                ENDDO

             ENDDO
          ENDIF

          DO l = 0, 3

             DO  m = 1, surf_lsm_v(l)%ns
!
!--             Spectral albedos for vegetation/pavement/water surfaces
                DO  ind_type = 0, 2
                   IF ( surf_lsm_v(l)%albedo_type(ind_type,m) /= 0 )  THEN
                      surf_lsm_v(l)%aldif(ind_type,m) =                        &
                            albedo_pars(0,surf_lsm_v(l)%albedo_type(ind_type,m))
                      surf_lsm_v(l)%asdif(ind_type,m) =                        &
                            albedo_pars(1,surf_lsm_v(l)%albedo_type(ind_type,m))
                      surf_lsm_v(l)%aldir(ind_type,m) =                        &
                            albedo_pars(0,surf_lsm_v(l)%albedo_type(ind_type,m))
                      surf_lsm_v(l)%asdir(ind_type,m) =                        &
                            albedo_pars(1,surf_lsm_v(l)%albedo_type(ind_type,m))
                      surf_lsm_v(l)%albedo(ind_type,m) =                       &
                            albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_type,m))
                   ENDIF
                ENDDO
             ENDDO
!
!--          For urban surface only if albedo has not been already initialized 
!--          in the urban-surface model via the ASCII file.
             IF ( .NOT. surf_usm_v(l)%albedo_from_ascii )  THEN
                DO  m = 1, surf_usm_v(l)%ns
!
!--                Spectral albedos for wall/green/window surfaces
                   DO  ind_type = 0, 2
                      IF ( surf_usm_v(l)%albedo_type(ind_type,m) /= 0 )  THEN
                         surf_usm_v(l)%aldif(ind_type,m) =                     &
                            albedo_pars(0,surf_usm_v(l)%albedo_type(ind_type,m))
                         surf_usm_v(l)%asdif(ind_type,m) =                     &
                            albedo_pars(1,surf_usm_v(l)%albedo_type(ind_type,m))
                         surf_usm_v(l)%aldir(ind_type,m) =                     &
                            albedo_pars(0,surf_usm_v(l)%albedo_type(ind_type,m))
                         surf_usm_v(l)%asdir(ind_type,m) =                     &
                            albedo_pars(1,surf_usm_v(l)%albedo_type(ind_type,m))
                         surf_usm_v(l)%albedo(ind_type,m) =                    &
                            albedo_pars(2,surf_usm_v(l)%albedo_type(ind_type,m))
                      ENDIF
                   ENDDO

                ENDDO
             ENDIF
          ENDDO
!
!--       Level 3 initialization at grid points where albedo type is zero.
!--       This case, spectral albedos are taken from file if available
          IF ( albedo_pars_f%from_file )  THEN
!
!--          Horizontal
             DO  m = 1, surf_lsm_h%ns
                i = surf_lsm_h%i(m)
                j = surf_lsm_h%j(m)
!
!--             Spectral albedos for vegetation/pavement/water surfaces
                DO  ind_type = 0, 2
                   IF ( surf_lsm_h%albedo_type(ind_type,m) == 0 )  THEN
                      IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )&
                         surf_lsm_h%albedo(ind_type,m) =                       &
                                                albedo_pars_f%pars_xy(1,j,i)
                      IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )&
                         surf_lsm_h%aldir(ind_type,m) =                        &
                                                albedo_pars_f%pars_xy(1,j,i)
                      IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )&
                         surf_lsm_h%aldif(ind_type,m) =                        &
                                                albedo_pars_f%pars_xy(2,j,i)
                      IF ( albedo_pars_f%pars_xy(3,j,i) /= albedo_pars_f%fill )&
                         surf_lsm_h%asdir(ind_type,m) =                        &
                                                albedo_pars_f%pars_xy(3,j,i)
                      IF ( albedo_pars_f%pars_xy(4,j,i) /= albedo_pars_f%fill )&
                         surf_lsm_h%asdif(ind_type,m) =                        &
                                                albedo_pars_f%pars_xy(4,j,i)
                   ENDIF
                ENDDO
             ENDDO
!
!--          For urban surface only if albedo has not been already initialized 
!--          in the urban-surface model via the ASCII file.
             IF ( .NOT. surf_usm_h%albedo_from_ascii )  THEN
                DO  m = 1, surf_usm_h%ns
                   i = surf_usm_h%i(m)
                   j = surf_usm_h%j(m)
!
!--                Spectral albedos for wall/green/window surfaces
                   DO  ind_type = 0, 2
                      IF ( surf_usm_h%albedo_type(ind_type,m) == 0 )  THEN
                         IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )&
                            surf_usm_h%albedo(ind_type,m) =                       &
                                                albedo_pars_f%pars_xy(1,j,i)
                         IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )&
                            surf_usm_h%aldir(ind_type,m) =                        &
                                                albedo_pars_f%pars_xy(1,j,i)
                         IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )&
                            surf_usm_h%aldif(ind_type,m) =                        &
                                                albedo_pars_f%pars_xy(2,j,i)
                         IF ( albedo_pars_f%pars_xy(3,j,i) /= albedo_pars_f%fill )&
                            surf_usm_h%asdir(ind_type,m) =                        &
                                                albedo_pars_f%pars_xy(3,j,i)
                         IF ( albedo_pars_f%pars_xy(4,j,i) /= albedo_pars_f%fill )&
                            surf_usm_h%asdif(ind_type,m) =                        &
                                                albedo_pars_f%pars_xy(4,j,i)
                      ENDIF
                   ENDDO

                ENDDO
             ENDIF
!
!--          Vertical
             DO  l = 0, 3
                ioff = surf_lsm_v(l)%ioff
                joff = surf_lsm_v(l)%joff

                DO  m = 1, surf_lsm_v(l)%ns
                   i = surf_lsm_v(l)%i(m)
                   j = surf_lsm_v(l)%j(m)
!
!--                Spectral albedos for vegetation/pavement/water surfaces
                   DO  ind_type = 0, 2
                      IF ( surf_lsm_v(l)%albedo_type(ind_type,m) == 0 )  THEN
                         IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /=        &
                              albedo_pars_f%fill )                             &
                            surf_lsm_v(l)%albedo(ind_type,m) =                 &
                                          albedo_pars_f%pars_xy(1,j+joff,i+ioff)
                         IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /=        &
                              albedo_pars_f%fill )                             &
                            surf_lsm_v(l)%aldir(ind_type,m) =                  &
                                          albedo_pars_f%pars_xy(1,j+joff,i+ioff)
                         IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /=        &
                              albedo_pars_f%fill )                             &
                            surf_lsm_v(l)%aldif(ind_type,m) =                  &
                                          albedo_pars_f%pars_xy(2,j+joff,i+ioff)
                         IF ( albedo_pars_f%pars_xy(3,j+joff,i+ioff) /=        &
                              albedo_pars_f%fill )                             &
                            surf_lsm_v(l)%asdir(ind_type,m) =                  &
                                          albedo_pars_f%pars_xy(3,j+joff,i+ioff)
                         IF ( albedo_pars_f%pars_xy(4,j+joff,i+ioff) /=        &
                              albedo_pars_f%fill )                             &
                            surf_lsm_v(l)%asdif(ind_type,m) =                  &
                                          albedo_pars_f%pars_xy(4,j+joff,i+ioff)
                      ENDIF
                   ENDDO
                ENDDO
!
!--             For urban surface only if albedo has not been already initialized 
!--             in the urban-surface model via the ASCII file.
                IF ( .NOT. surf_usm_v(l)%albedo_from_ascii )  THEN
                   ioff = surf_usm_v(l)%ioff
                   joff = surf_usm_v(l)%joff

                   DO  m = 1, surf_usm_v(l)%ns
                      i = surf_usm_v(l)%i(m)
                      j = surf_usm_v(l)%j(m)
!
!--                   Spectral albedos for wall/green/window surfaces
                      DO  ind_type = 0, 2
                         IF ( surf_usm_v(l)%albedo_type(ind_type,m) == 0 )  THEN
                            IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /=        &
                                 albedo_pars_f%fill )                             &
                               surf_usm_v(l)%albedo(ind_type,m) =                 &
                                             albedo_pars_f%pars_xy(1,j+joff,i+ioff)
                            IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /=        &
                                 albedo_pars_f%fill )                             &
                               surf_usm_v(l)%aldir(ind_type,m) =                  &
                                             albedo_pars_f%pars_xy(1,j+joff,i+ioff)
                            IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /=        &
                                 albedo_pars_f%fill )                             &
                               surf_usm_v(l)%aldif(ind_type,m) =                  &
                                             albedo_pars_f%pars_xy(2,j+joff,i+ioff)
                            IF ( albedo_pars_f%pars_xy(3,j+joff,i+ioff) /=        &
                                 albedo_pars_f%fill )                             &
                               surf_usm_v(l)%asdir(ind_type,m) =                  &
                                             albedo_pars_f%pars_xy(3,j+joff,i+ioff)
                            IF ( albedo_pars_f%pars_xy(4,j+joff,i+ioff) /=        &
                                 albedo_pars_f%fill )                             &
                               surf_usm_v(l)%asdif(ind_type,m) =                  &
                                             albedo_pars_f%pars_xy(4,j+joff,i+ioff)
                         ENDIF
                      ENDDO

                   ENDDO
                ENDIF
             ENDDO

          ENDIF

!
!--       Calculate initial values of current (cosine of) the zenith angle and 
!--       whether the sun is up
          CALL calc_zenith
!
!--       readjust date and time to its initial value
          CALL init_date_and_time
!
!--       Calculate initial surface albedo for different surfaces
          IF ( .NOT. constant_albedo )  THEN
#if defined( __netcdf )
!
!--          Horizontally aligned natural and urban surfaces
             CALL calc_albedo( surf_lsm_h )
             CALL calc_albedo( surf_usm_h )
!
!--          Vertically aligned natural and urban surfaces
             DO  l = 0, 3
                CALL calc_albedo( surf_lsm_v(l) )
                CALL calc_albedo( surf_usm_v(l) )
             ENDDO
#endif
          ELSE
!
!--          Initialize sun-inclination independent spectral albedos
!--          Horizontal surfaces
             IF ( surf_lsm_h%ns > 0 )  THEN
                surf_lsm_h%rrtm_aldir = surf_lsm_h%aldir
                surf_lsm_h%rrtm_asdir = surf_lsm_h%asdir
                surf_lsm_h%rrtm_aldif = surf_lsm_h%aldif
                surf_lsm_h%rrtm_asdif = surf_lsm_h%asdif
             ENDIF
             IF ( surf_usm_h%ns > 0 )  THEN
                surf_usm_h%rrtm_aldir = surf_usm_h%aldir
                surf_usm_h%rrtm_asdir = surf_usm_h%asdir
                surf_usm_h%rrtm_aldif = surf_usm_h%aldif
                surf_usm_h%rrtm_asdif = surf_usm_h%asdif
             ENDIF
!
!--          Vertical surfaces
             DO  l = 0, 3
                IF ( surf_lsm_v(l)%ns > 0 )  THEN
                   surf_lsm_v(l)%rrtm_aldir = surf_lsm_v(l)%aldir
                   surf_lsm_v(l)%rrtm_asdir = surf_lsm_v(l)%asdir
                   surf_lsm_v(l)%rrtm_aldif = surf_lsm_v(l)%aldif
                   surf_lsm_v(l)%rrtm_asdif = surf_lsm_v(l)%asdif
                ENDIF
                IF ( surf_usm_v(l)%ns > 0 )  THEN
                   surf_usm_v(l)%rrtm_aldir = surf_usm_v(l)%aldir
                   surf_usm_v(l)%rrtm_asdir = surf_usm_v(l)%asdir
                   surf_usm_v(l)%rrtm_aldif = surf_usm_v(l)%aldif
                   surf_usm_v(l)%rrtm_asdif = surf_usm_v(l)%asdif
                ENDIF
             ENDDO

          ENDIF

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

!
!--       Read trace gas profiles from file. This routine provides
!--       the rrtm_ arrays (1:nzt_rad+1)
          CALL read_trace_gas_data
#endif
       ENDIF

!
!--    Perform user actions if required
       CALL user_init_radiation

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

          CASE ( 'rrtmg' )
             CALL radiation_rrtmg

          CASE ( 'clear-sky' )
             CALL radiation_clearsky

          CASE ( 'constant' )
             CALL radiation_constant

          CASE DEFAULT

       END SELECT

! readjust date and time to its initial value
       CALL init_date_and_time

       CALL location_message( 'finished', .TRUE. )

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

!
!--    If required, read or calculate and write out the SVF 
       IF ( radiation_interactions .AND. read_svf)  THEN
!
!--       Read sky-view factors and further required data from file
          CALL location_message( '    Start reading SVF from file', .FALSE. )
          CALL radiation_read_svf()
          CALL location_message( '    Reading SVF from file has finished', .TRUE. )

       ELSEIF ( radiation_interactions .AND. .NOT. read_svf)  THEN
!
!--       calculate SFV and CSF
          CALL location_message( '    Start calculation of SVF', .FALSE. )
          CALL radiation_calc_svf()
          CALL location_message( '    Calculation of SVF has finished', .TRUE. )
       ENDIF

       IF ( radiation_interactions .AND. write_svf)  THEN
!
!--       Write svf, csf svfsurf and csfsurf data to file
          CALL location_message( '    Start writing SVF in file', .FALSE. )
          CALL radiation_write_svf()
          CALL location_message( '    Writing SVF in file has finished', .TRUE. )
       ENDIF

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

       RETURN

    END SUBROUTINE radiation_init


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


       IMPLICIT NONE

       INTEGER(iwp) ::  l         !< running index for surface orientation
       REAL(wp)     ::  pt1       !< potential temperature at first grid level or mean value at urban layer top 
       REAL(wp)     ::  pt1_l     !< potential temperature at first grid level or mean value at urban layer top at local subdomain 
       REAL(wp)     ::  ql1       !< liquid water mixing ratio at first grid level or mean value at urban layer top 
       REAL(wp)     ::  ql1_l     !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain 

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

!
!--    Calculate current zenith angle
       CALL calc_zenith

!
!--    Calculate sky transmissivity
       sky_trans = 0.6_wp + 0.2_wp * cos_zenith

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

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

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

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

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

       CONTAINS

          SUBROUTINE radiation_clearsky_surf

             IMPLICIT NONE

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

             IF ( surf%ns < 1 )  RETURN

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

                k = nz_urban_t

                surf%rad_sw_in  = solar_constant * sky_trans * cos_zenith
                surf%rad_sw_out = albedo_urb * surf%rad_sw_in
                
                surf%rad_lw_in  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k+1))**4

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

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

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

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

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

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

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

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

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


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

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

                ENDDO

             ENDIF

!
!--          Fill out values in radiation arrays
             DO  m = 1, surf%ns
                i = surf%i(m)
                j = surf%j(m)
                rad_sw_in(0,j,i) = surf%rad_sw_in(m)
                rad_sw_out(0,j,i) = surf%rad_sw_out(m)
                rad_lw_in(0,j,i) = surf%rad_lw_in(m)
                rad_lw_out(0,j,i) = surf%rad_lw_out(m)
             ENDDO
 
          END SUBROUTINE radiation_clearsky_surf

    END SUBROUTINE radiation_clearsky


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


       IMPLICIT NONE

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

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

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

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

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

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

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

       CONTAINS

          SUBROUTINE radiation_constant_surf

             IMPLICIT NONE

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

             IF ( surf%ns < 1 )  RETURN

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

                surf%rad_net = net_radiation

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

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

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

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

                surf%rad_sw_out =  albedo_urb * surf%rad_sw_in

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

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

                   surf%rad_net(m) = net_radiation

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

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

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

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

                ENDDO

             ENDIF

!
!--          Fill out values in radiation arrays
             DO  m = 1, surf%ns
                i = surf%i(m)
                j = surf%j(m)
                rad_sw_in(0,j,i) = surf%rad_sw_in(m)
                rad_sw_out(0,j,i) = surf%rad_sw_out(m)
                rad_lw_in(0,j,i) = surf%rad_lw_in(m)
                rad_lw_out(0,j,i) = surf%rad_lw_out(m)
             ENDDO

          END SUBROUTINE radiation_constant_surf
          

    END SUBROUTINE radiation_constant

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


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

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

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

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

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


    END SUBROUTINE radiation_header
   

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


       IMPLICIT NONE

       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file 
       
       NAMELIST /radiation_par/   albedo, albedo_lw_dif, albedo_lw_dir,         &
                                  albedo_sw_dif, albedo_sw_dir, albedo_type,    &
                                  constant_albedo, dt_radiation, emissivity,    &
                                  lw_radiation, max_raytracing_dist,            &
                                  min_irrf_value, mrt_geom_human,               &
                                  mrt_include_sw, mrt_nlevels,                  &
                                  mrt_skip_roof, net_radiation, nrefsteps,      &
                                  plant_lw_interact, rad_angular_discretization,&
                                  radiation_interactions_on, radiation_scheme,  &
                                  raytrace_discrete_azims,                      &
                                  raytrace_discrete_elevs, raytrace_mpi_rma,    &
                                  skip_time_do_radiation, surface_reflections,  &
                                  svfnorm_report_thresh, sw_radiation,          &
                                  unscheduled_radiation_calls

   
       NAMELIST /radiation_parameters/ albedo, albedo_lw_dif, albedo_lw_dir,    &
                                  albedo_sw_dif, albedo_sw_dir, albedo_type,    &
                                  constant_albedo, dt_radiation, emissivity,    &
                                  lw_radiation, max_raytracing_dist,            &
                                  min_irrf_value, mrt_geom_human,               &
                                  mrt_include_sw, mrt_nlevels,                  &
                                  mrt_skip_roof, net_radiation, nrefsteps,      &
                                  plant_lw_interact, rad_angular_discretization,&
                                  radiation_interactions_on, radiation_scheme,  &
                                  raytrace_discrete_azims,                      &
                                  raytrace_discrete_elevs, raytrace_mpi_rma,    &
                                  skip_time_do_radiation, surface_reflections,  &
                                  svfnorm_report_thresh, sw_radiation,          &
                                  unscheduled_radiation_calls
   
       line = ' '
       
!
!--    Try to find radiation model namelist
       REWIND ( 11 )
       line = ' '
       DO WHILE ( INDEX( line, '&radiation_parameters' ) == 0 )
          READ ( 11, '(A)', END=12 )  line
       ENDDO
       BACKSPACE ( 11 )

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

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

       GOTO 14

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

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

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

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

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

       GOTO 14

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

 14    CONTINUE
       
    END SUBROUTINE radiation_parin


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

#if defined ( __rrtmg )
       USE indices,                                                            &
           ONLY:  nbgp

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

       IMPLICIT NONE


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

       REAL(wp)     ::  nc_rad, &    !< number concentration of cloud droplets
                        s_r2,   &    !< weighted sum over all droplets with r^2
                        s_r3         !< weighted sum over all droplets with r^3

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

!
!--    Calculate current (cosine of) zenith angle and whether the sun is up
       CALL calc_zenith     
       zenith(0) = cos_zenith
!
!--    Calculate surface albedo. In case average radiation is applied, 
!--    this is not required.
#if defined( __netcdf )
       IF ( .NOT. constant_albedo )  THEN
!
!--       Horizontally aligned default, natural and urban surfaces
          CALL calc_albedo( surf_lsm_h    )
          CALL calc_albedo( surf_usm_h    )
!
!--       Vertically aligned default, natural and urban surfaces
          DO  l = 0, 3
             CALL calc_albedo( surf_lsm_v(l) )
             CALL calc_albedo( surf_usm_v(l) )
          ENDDO
       ENDIF
#endif

!
!--    Prepare input data for RRTMG

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


       IF ( average_radiation ) THEN

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

          rrtm_emis = emissivity_urb
!
!--       Calculate mean pt profile. Actually, only one height level is required.
          CALL calc_mean_profile( pt, 4 )
          pt_av = hom(:, 1, 4, 0)
          
          IF ( humidity )  THEN
             CALL calc_mean_profile( q, 41 )
             q_av  = hom(:, 1, 41, 0)
          ENDIF
!
!--       Prepare profiles of temperature and H2O volume mixing ratio
          rrtm_tlev(0,nzb+1) = t_rad_urb

          IF ( bulk_cloud_model )  THEN

             CALL calc_mean_profile( ql, 54 )
             ! average ql is now in hom(:, 1, 54, 0)
             ql_av = hom(:, 1, 54, 0)
             
             DO k = nzb+1, nzt+1
                rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp       &
                                 )**.286_wp + lv_d_cp * ql_av(k)
                rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q_av(k) - ql_av(k))
             ENDDO
          ELSE
             DO k = nzb+1, nzt+1
                rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp       &
                                 )**.286_wp
             ENDDO

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

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

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

          ENDDO

!--       Linear interpolate to zw grid. Loop reaches one level further up
!--       due to the staggered grid in RRTMG
          DO k = nzb+2, nzt+8
             rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) -        &
                                rrtm_tlay(0,k-1))                           &
                                / ( rrtm_play(0,k) - rrtm_play(0,k-1) )     &
                                * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
          ENDDO


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

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

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

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

!
!--       Set surface temperature
          rrtm_tsfc = t_rad_urb
          
          IF ( lw_radiation )  THEN        
          
             CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
             rrtm_play       , rrtm_plev    , rrtm_tlay    , rrtm_tlev      ,&
             rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr   , rrtm_co2vmr    ,&
             rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr   , rrtm_cfc11vmr  ,&
             rrtm_cfc12vmr   , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis      ,&
             rrtm_inflglw    , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr     ,&
             rrtm_lw_taucld  , rrtm_cicewp  , rrtm_cliqwp  , rrtm_reice     ,& 
             rrtm_reliq      , rrtm_lw_tauaer,                               &
             rrtm_lwuflx     , rrtm_lwdflx  , rrtm_lwhr  ,                   &
             rrtm_lwuflxc    , rrtm_lwdflxc , rrtm_lwhrc ,                   &
             rrtm_lwuflx_dt  ,  rrtm_lwuflxc_dt )

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

          ENDIF

          IF ( sw_radiation .AND. sun_up )  THEN
             CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld     , rrtm_iaer      ,&
             rrtm_play      , rrtm_plev     , rrtm_tlay     , rrtm_tlev      ,&
             rrtm_tsfc      , rrtm_h2ovmr   , rrtm_o3vmr    , rrtm_co2vmr    ,&
             rrtm_ch4vmr    , rrtm_n2ovmr   , rrtm_o2vmr    , rrtm_asdir     ,&
             rrtm_asdif     , rrtm_aldir    , rrtm_aldif    , zenith         ,&
             0.0_wp         , day_of_year   , solar_constant, rrtm_inflgsw   ,&
             rrtm_iceflgsw  , rrtm_liqflgsw , rrtm_cldfr    , rrtm_sw_taucld ,&
             rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp    ,&
             rrtm_cliqwp    , rrtm_reice    , rrtm_reliq    , rrtm_sw_tauaer ,&
             rrtm_sw_ssaaer , rrtm_sw_asmaer, rrtm_sw_ecaer , rrtm_swuflx    ,&
             rrtm_swdflx    , rrtm_swhr     , rrtm_swuflxc  , rrtm_swdflxc   ,&
             rrtm_swhrc     , rrtm_dirdflux , rrtm_difdflux )
  
!
!--          Save fluxes:
!--          - whole domain
             DO k = nzb, nzt+1
                rad_sw_in(k,:,:)  = rrtm_swdflx(0,k)
                rad_sw_out(k,:,:) = rrtm_swuflx(0,k)
             ENDDO
!--          - direct and diffuse SW at urban-surface-layer (required by RTM)
             rad_sw_in_dir(:,:) = rrtm_dirdflux(0,nzb)
             rad_sw_in_diff(:,:) = rrtm_difdflux(0,nzb)

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

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


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

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

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

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

                ENDDO

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


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

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

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

!
!--                      Calculate cloud droplet effective radius
                         IF ( bulk_cloud_model )  THEN
!
!--                         Calculete effective droplet radius. In case of using 
!--                         cloud_scheme = 'morrison' and a non reasonable number 
!--                         of cloud droplets the inital aerosol number  
!--                         concentration is considered.
                            IF ( microphysics_morrison )  THEN
                               IF ( nc(k,j,i) > 1.0E-20_wp )  THEN
                                  nc_rad = nc(k,j,i)
                               ELSE
                                  nc_rad = na_init
                               ENDIF 
                            ELSE
                               nc_rad = nc_const
                            ENDIF  

                            rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)     &
                                              * rho_surface                       &
                                              / ( 4.0_wp * pi * nc_rad * rho_l )  &
                                              )**0.33333333333333_wp              &
                                              * EXP( LOG( sigma_gc )**2 )

                         ELSEIF ( cloud_droplets )  THEN
                            number_of_particles = prt_count(k,j,i)

                            IF (number_of_particles <= 0)  CYCLE
                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
                            s_r2 = 0.0_wp
                            s_r3 = 0.0_wp

                            DO  n = 1, number_of_particles
                               IF ( particles(n)%particle_mask )  THEN
                                  s_r2 = s_r2 + particles(n)%radius**2 *       &
                                         particles(n)%weight_factor
                                  s_r3 = s_r3 + particles(n)%radius**3 *       &
                                         particles(n)%weight_factor
                               ENDIF
                            ENDDO

                            IF ( s_r2 > 0.0_wp )  rrtm_reliq(0,k) = s_r3 / s_r2

                         ENDIF

!
!--                      Limit effective radius
                         IF ( rrtm_reliq(0,k) > 0.0_wp )  THEN
                            rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
                            rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
                        ENDIF
                      ENDIF
                   ENDDO
                ENDIF

!
!--             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.  
!--             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. 
                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
                   rrtm_emis = surf_lsm_h%frac(ind_veg_wall,m)  *              &
                               surf_lsm_h%emissivity(ind_veg_wall,m)  +        &
                               surf_lsm_h%frac(ind_pav_green,m) *              &
                               surf_lsm_h%emissivity(ind_pav_green,m) +        & 
                               surf_lsm_h%frac(ind_wat_win,m)   *              &
                               surf_lsm_h%emissivity(ind_wat_win,m)
                   rrtm_tsfc = surf_lsm_h%pt_surface(m) * exner(nzb)
                ENDDO             
                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
                   rrtm_emis = surf_usm_h%frac(ind_veg_wall,m)  *              &
                               surf_usm_h%emissivity(ind_veg_wall,m)  +        &
                               surf_usm_h%frac(ind_pav_green,m) *              &
                               surf_usm_h%emissivity(ind_pav_green,m) +        & 
                               surf_usm_h%frac(ind_wat_win,m)   *              &
                               surf_usm_h%emissivity(ind_wat_win,m)
                   rrtm_tsfc = surf_usm_h%pt_surface(m) * exner(nzb)
                ENDDO
!
!--             Obtain topography top index (lower bound of RRTMG)
                k_topo = get_topography_top_index_ji( j, i, 's' )

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

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

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

                   DEALLOCATE ( rrtm_lw_taucld_dum )
                   DEALLOCATE ( rrtm_lw_tauaer_dum )
!
!--                Save fluxes
                   DO k = k_topo, nzt+1
                      rad_lw_in(k,j,i)  = rrtm_lwdflx(0,k)
                      rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
                   ENDDO

!
!--                Save heating rates (convert from K/d to K/h)
                   DO k = k_topo+1, nzt+1
                      rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k-k_topo)  * d_hours_day
                      rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k-k_topo) * d_hours_day
                   ENDDO

!
!--                Save surface radiative fluxes and change in LW heating rate 
!--                onto respective surface elements
!--                Horizontal surfaces
                   DO  m = surf_lsm_h%start_index(j,i),                        &
                           surf_lsm_h%end_index(j,i)
                      surf_lsm_h%rad_lw_in(m)           = rrtm_lwdflx(0,k_topo)
                      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             
                   DO  m = surf_usm_h%start_index(j,i),                        &
                           surf_usm_h%end_index(j,i)
                      surf_usm_h%rad_lw_in(m)           = rrtm_lwdflx(0,k_topo)
                      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 
!--                respective surface element
                   DO  l = 0, 3
                      DO  m = surf_lsm_v(l)%start_index(j,i),                  &
                              surf_lsm_v(l)%end_index(j,i)
                         k                                    = surf_lsm_v(l)%k(m)
                         surf_lsm_v(l)%rad_lw_in(m)           = rrtm_lwdflx(0,k)
                         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             
                      DO  m = surf_usm_v(l)%start_index(j,i),                  &
                              surf_usm_v(l)%end_index(j,i)
                         k                                    = surf_usm_v(l)%k(m)
                         surf_usm_v(l)%rad_lw_in(m)           = rrtm_lwdflx(0,k)
                         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

                ENDIF

                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 
!--                column model
!--                (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) )
                      rrtm_asdif(1)  = SUM( surf_lsm_h%frac(:,m) *             &
                                            surf_lsm_h%rrtm_asdif(:,m) )
                      rrtm_aldir(1)  = SUM( surf_lsm_h%frac(:,m) *             &
                                            surf_lsm_h%rrtm_aldir(:,m) )
                      rrtm_aldif(1)  = SUM( surf_lsm_h%frac(:,m) *             &
                                            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) )
                      rrtm_asdif(1)  = SUM( surf_usm_h%frac(:,m) *             &
                                            surf_usm_h%rrtm_asdif(:,m) )
                      rrtm_aldir(1)  = SUM( surf_usm_h%frac(:,m) *             &
                                            surf_usm_h%rrtm_aldir(:,m) )
                      rrtm_aldif(1)  = SUM( surf_usm_h%frac(:,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 
!--                rrtmg_sw is called.
!--                As one dimesion is allocated with zero size, compiler complains
!--                that rank of the array does not match that of the 
!--                assumed-shaped arguments in the RRTMG library. In order to 
!--                avoid this, write to dummy arguments and give pass the entire 
!--                dummy array. Seems to be the only existing work-around.  
                   ALLOCATE( rrtm_sw_taucld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
                   ALLOCATE( rrtm_sw_ssacld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
                   ALLOCATE( rrtm_sw_asmcld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
                   ALLOCATE( rrtm_sw_fsfcld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
                   ALLOCATE( rrtm_sw_tauaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
                   ALLOCATE( rrtm_sw_ssaaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
                   ALLOCATE( rrtm_sw_asmaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
                   ALLOCATE( rrtm_sw_ecaer_dum(0:0,k_topo+1:nzt_rad+1,1:naerec+1)  )
     
                   rrtm_sw_taucld_dum = rrtm_sw_taucld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
                   rrtm_sw_ssacld_dum = rrtm_sw_ssacld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
                   rrtm_sw_asmcld_dum = rrtm_sw_asmcld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
                   rrtm_sw_fsfcld_dum = rrtm_sw_fsfcld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
                   rrtm_sw_tauaer_dum = rrtm_sw_tauaer(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1)
                   rrtm_sw_ssaaer_dum = rrtm_sw_ssaaer(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1)
                   rrtm_sw_asmaer_dum = rrtm_sw_asmaer(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1)
                   rrtm_sw_ecaer_dum  = rrtm_sw_ecaer(0:0,k_topo+1:nzt_rad+1,1:naerec+1)

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

                   DEALLOCATE( rrtm_sw_taucld_dum )
                   DEALLOCATE( rrtm_sw_ssacld_dum )
                   DEALLOCATE( rrtm_sw_asmcld_dum )
                   DEALLOCATE( rrtm_sw_fsfcld_dum )
                   DEALLOCATE( rrtm_sw_tauaer_dum )
                   DEALLOCATE( rrtm_sw_ssaaer_dum )
                   DEALLOCATE( rrtm_sw_asmaer_dum )
                   DEALLOCATE( rrtm_sw_ecaer_dum )
!
!--                Save fluxes
                   DO k = nzb, nzt+1
                      rad_sw_in(k,j,i)  = rrtm_swdflx(0,k)
                      rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
                   ENDDO
!
!--                Save heating rates (convert from K/d to K/s)
                   DO k = nzb+1, nzt+1
                      rad_sw_hr(k,j,i)     = rrtm_swhr(0,k)  * d_hours_day
                      rad_sw_cs_hr(k,j,i)  = rrtm_swhrc(0,k) * d_hours_day
                   ENDDO

!
!--                Save surface radiative fluxes onto respective surface elements
!--                Horizontal surfaces
                   DO  m = surf_lsm_h%start_index(j,i),                        &
                           surf_lsm_h%end_index(j,i)
                      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             
                   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 
!--                level of the surface element
                   DO  l = 0, 3
                      DO  m = surf_lsm_v(l)%start_index(j,i),                  &
                              surf_lsm_v