Index: palm/trunk/SOURCE/radiation_model.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model.f90	(revision 1572)
+++ palm/trunk/SOURCE/radiation_model.f90	(revision 1585)
@@ -15,10 +15,10 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Added support for RRTMG
 ! 
 ! Former revisions:
@@ -40,68 +40,262 @@
 ! Description:
 ! ------------
-! Radiation model(s), to be used e.g. with the land surface scheme
+! Radiation models and interfaces
+! To do: move variable definitions used in init_radiation only to the subroutine
+! as they are no longer required after initialization.
+! Note that many variables have a leading dummy dimension (0:0) in order to
+! match the assume-size shape expected by the RRTMG model
+!
+! To do:
+! - Output of full column vertical profiles used in RRTMG
+! - Output of other rrtm arrays (such as volume mixing ratios)
+! - Adapt for use with topography
+! - Remove 3D dummy arrays (such as clear-sky output)
 !------------------------------------------------------------------------------!
 
     USE arrays_3d,                                                             &
-        ONLY: pt
+        ONLY:  hyp, pt, q, ql, zw
+
+    USE cloud_parameters,                                                      &
+        ONLY:  cp, l_v, nc_const, rho_l, sigma_gc  
+
+    USE constants,                                                             &
+        ONLY:  pi
 
     USE control_parameters,                                                    &
-        ONLY: phi, surface_pressure, time_since_reference_point
+        ONLY:  cloud_droplets, cloud_physics, g, initializing_actions,         &
+               large_scale_forcing, lsf_surf, phi, pt_surface,                 &
+               surface_pressure, time_since_reference_point
 
     USE indices,                                                               &
-        ONLY:  nxlg, nxrg, nyng, nysg, nzb_s_inner
+        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb_s_inner, nzb, nzt
 
     USE kinds
+
+    USE netcdf
 
     USE netcdf_control,                                                        &
         ONLY:  dots_label, dots_num, dots_unit
 
+#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
+
+
 
     IMPLICIT NONE
 
-    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtm'
-
-    INTEGER(iwp) :: i, j, k
-
-
-    INTEGER(iwp) :: day_init     = 172,  & !: day of the year at model start (21/06)
-                    dots_rad     = 0,    & !: starting index for timeseries output 
-                    irad_scheme  = 0
-
-    LOGICAL :: radiation = .FALSE.  !: flag parameter indicating wheather the radiation model is used
-
-    REAL(wp), PARAMETER :: sw_0 = 1368.0_wp, &    !: solar constant  
-                           pi = 3.14159265358979323_wp, &
-                           sigma_sb  = 5.67E-8_wp !: Stefan-Boltzmann constant
+    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
+
+!
+!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
+    CHARACTER(37), DIMENSION(0:16), 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
+                                   'landice',                              & ! 14 
+                                   'sea ice',                              & ! 15 
+                                   'snow'                                  & ! 16
+                                                         /)
+
+    INTEGER(iwp) :: albedo_type  = 5,    & !: Albedo surface type (default: short grassland)
+                    day,                 & !: current day of the year
+                    day_init     = 172,  & !: day of the year at model start (21/06)
+                    dots_rad     = 0       !: starting index for timeseries output 
+
+
+
+
+
+
+    LOGICAL ::  constant_albedo = .FALSE., & !: flag parameter indicating whether the albedo may change depending on zenith
+                lw_radiation = .TRUE.,     & !: flag parameter indicing 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 indicing whether shortwave radiation shall be calculated
+
+    REAL(wp), PARAMETER :: sigma_sb       = 5.67E-8_wp, & !: Stefan-Boltzmann constant
+                           solar_constant = 1368.0_wp     !: solar constant at top of atmosphere
  
-    REAL(wp) :: albedo = 0.2_wp,             & !: NAMELIST alpha
-                dt_radiation = 0.0_wp,       & !: radiation model timestep
-                exn,                         & !: Exner function
-                lon = 0.0_wp,                & !: longitude in radians
-                lat = 0.0_wp,                & !: latitude in radians
-                decl_1,                      & !: declination coef. 1
-                decl_2,                      & !: declination coef. 2
-                decl_3,                      & !: declination coef. 3
-                time_utc,                    & !: current time in UTC
-                time_utc_init = 43200.0_wp,  & !: UTC time at model start (noon)
-                day,                         & !: current day of the year
-                lambda = 0.0_wp,             & !: longitude in degrees
-                declination,                 & !: solar declination angle
-                net_radiation = 0.0_wp,      & !: net radiation at surface
-                hour_angle,                  & !: solar hour angle
-                time_radiation = 0.0_wp,     &
-                zenith,                      & !: solar zenith angle
-                sky_trans                      !: sky transmissivity
+    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
+                dt_radiation = 0.0_wp,        & !: radiation model timestep
+                emissivity = 0.95_wp,         & !: NAMELIST surface emissivity
+                exn,                          & !: Exner function
+                lon = 0.0_wp,                 & !: longitude in radians
+                lat = 0.0_wp,                 & !: latitude in radians
+                decl_1,                       & !: declination coef. 1
+                decl_2,                       & !: declination coef. 2
+                decl_3,                       & !: declination coef. 3
+                time_utc,                     & !: current time in UTC
+                time_utc_init = 43200.0_wp,   & !: UTC time at model start (noon)
+                lambda = 0.0_wp,              & !: longitude in degrees
+                net_radiation = 0.0_wp,       & !: net radiation at surface
+                time_radiation = 0.0_wp,      & !: time since last call of radiation code
+                sky_trans                       !: sky transmissivity
+
+
+    REAL(wp), DIMENSION(0:0) ::  zenith        !: solar zenith angle
 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
-                alpha,                       & !: surface albedo
+                alpha,                       & !: surface broadband albedo (used for clear-sky scheme)
                 rad_net,                     & !: net radiation at the surface
-                rad_net_av,                  & !: average of rad_net
-                rad_lw_in,                   & !: incoming longwave radiation
-                rad_lw_out,                  & !: outgoing longwave radiation
-                rad_sw_in,                   & !: incoming shortwave radiation
-                rad_sw_in_av,                & !: average of rad_sw_in
-                rad_sw_out                     !: outgoing shortwave radiation
-
+                rad_net_av                     !: average of rad_net
+
+!
+!-- Land surface albedos for solar zenith angle of 60Â° after Briegleb (1992)     
+!-- (shortwave, longwave, broadband):   sw,      lw,      bb,
+    REAL(wp), DIMENSION(0:2,1:15), 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.95_wp, 0.70_wp, 0.82_wp             & ! 15
+                                 /), (/ 3, 15 /) )
+
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
+                        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
+                        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
+
+!
+!-- 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 for RRTMGS (should not be changed)
+    INTEGER(iwp), PARAMETER :: 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)
+                    rrtm_idrv = 0        !: longwave upward flux calculation option (0,1)
+
+    LOGICAL :: snd_exists = .FALSE.      !: flag parameter to check whether a user-defined input files exists
+
+    REAL(wp), PARAMETER :: mol_weight_air_d_wv = 1.607793_wp !: molecular weight dry air / water vapor
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !: hypostatic pressure from sounding data (hPa)
+                                           q_snd,       & !: specific humidity from sounding data (kg/kg) - dummy at the moment
+                                           rrtm_tsfc,   & !: dummy array for storing surface temperature
+                                           t_snd          !: actual temperature from sounding data (hPa)
+
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif,         & !: longwave diffuse albedo solar angle of 60Â°
+                                             aldir,         & !: longwave direct albedo solar angle of 60Â°
+                                             asdif,         & !: shortwave diffuse albedo solar angle of 60Â°
+                                             asdir,         & !: shortwave direct albedo solar angle of 60Â°
+                                             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/mÂ²)
+                                             rrtm_cldfr,    & !: cloud fraction (0,1)
+                                             rrtm_cliqwp,   & !: in-cloud liquid water path (g/mÂ²)
+                                             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_lwuflx,   & !: RRTM output of outgoing longwave radiation flux (W/m2)
+                                             rrtm_lwhr,     & !: RRTM output of longwave radiation heating rate (K/d)
+                                             rrtm_lwuflxc,  & !: RRTM output of incoming clear sky longwave radiation flux (W/m2)
+                                             rrtm_lwdflxc,  & !: RRTM output of outgoing clear sky longwave radiation flux (W/m2) 
+                                             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_swuflx,   & !: RRTM output of outgoing shortwave radiation flux (W/m2)
+                                             rrtm_swhr,     & !: RRTM output of shortwave radiation heating rate (K/d)
+                                             rrtm_swuflxc,  & !: RRTM output of incoming clear sky shortwave radiation flux (W/m2)
+                                             rrtm_swdflxc,  & !: RRTM output of outgoing clear sky shortwave radiation flux (W/m2) 
+                                             rrtm_swhrc       !: RRTM output of incoming shortwave clear sky radiation heating rate (K/d) 
+
+!
+!-- 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_lw_cs_hr,   & !: longwave clear sky radiation heating rate (K/d) (not used)
+                                                rad_lw_hr,      & !: longwave radiation heating rate (K/d)
+                                                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)
+                                                rad_sw_cs_hr,   & !: shortwave clear sky radiation heating rate (K/d) (not used)
+                                                rad_sw_hr,      & !: shortwave radiation heating rate (K/d)
+                                                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
+                                                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
 
     INTERFACE init_radiation
@@ -113,12 +307,12 @@
     END INTERFACE radiation_clearsky
 
-    INTERFACE radiation_constant
-       MODULE PROCEDURE radiation_constant
-    END INTERFACE radiation_constant
-
-    INTERFACE radiation_rrtm
-       MODULE PROCEDURE radiation_rrtm
-    END INTERFACE radiation_rrtm
-
+    INTERFACE radiation_rrtmg
+       MODULE PROCEDURE radiation_rrtmg
+    END INTERFACE radiation_rrtmg
+
+    INTERFACE radiation_tendency
+       MODULE PROCEDURE radiation_tendency
+       MODULE PROCEDURE radiation_tendency_ij
+    END INTERFACE radiation_tendency
 
     SAVE
@@ -126,11 +320,16 @@
     PRIVATE
 
-    PUBLIC albedo, day_init, dots_rad, dt_radiation, init_radiation,           &
-           irad_scheme, lambda, net_radiation, rad_net, rad_net_av, radiation, &
-           radiation_clearsky, radiation_constant, radiation_rrtm,             &
-           radiation_scheme, rad_sw_in, rad_sw_in_av, sigma_sb,                &
-           time_radiation, time_utc_init 
-
-
+    PUBLIC albedo, albedo_type, albedo_type_name, albedo_lw_dif, albedo_lw_dir,&
+           albedo_sw_dif, albedo_sw_dir, constant_albedo, day_init, dots_rad,  &
+           dt_radiation, init_radiation, lambda, lw_radiation, net_radiation,  &
+           rad_net, rad_net_av, radiation, radiation_clearsky,                 &
+           radiation_rrtmg, radiation_scheme, radiation_tendency, rad_lw_in,   &
+           rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_sw_in, rad_sw_in_av,   &
+           rad_sw_out, rad_sw_out_av, sigma_sb, sw_radiation, time_radiation,  &
+           time_utc_init 
+
+#if defined ( __rrtmg )
+    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
+#endif
 
  CONTAINS
@@ -143,45 +342,249 @@
     SUBROUTINE init_radiation
     
-
        IMPLICIT NONE
 
-       ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( rad_lw_in(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( rad_lw_out(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( rad_sw_in(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( rad_sw_out(nysg:nyng,nxlg:nxrg) )
-
-       rad_sw_in  = 0.0_wp
-       rad_sw_out = 0.0_wp
-       rad_lw_in  = 0.0_wp
-       rad_lw_out = 0.0_wp
-       rad_net    = 0.0_wp
-
-       alpha = albedo
+!
+!--    Allocate array for storing the surface net radiation
+       IF ( .NOT. ALLOCATED ( rad_net ) )  THEN
+          ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
+          rad_net = 0.0_wp
+       ENDIF
 
 !
 !--    Fix net radiation in case of radiation_scheme = 'constant'
-       IF ( irad_scheme == 0 )  THEN
+       IF ( radiation_scheme == 'constant' )  THEN
           rad_net = net_radiation
-!
-!--    Calculate radiation scheme constants
-       ELSEIF ( irad_scheme == 1 .OR. irad_scheme == 2 )  THEN
+          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
-!
-!--       Calculate latitude and longitude angles (lat and lon, respectively)
-          lat = phi * pi / 180.0_wp
-          lon = lambda * pi / 180.0_wp
+          lat    = phi * pi / 180.0_wp
+          lon    = lambda * pi / 180.0_wp
        ENDIF
+
+
+       IF ( radiation_scheme == 'clear-sky' )  THEN
+
+          ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
+
+          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_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 ) )  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
+
+          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
+
+          rad_sw_in  = 0.0_wp
+          rad_sw_out = 0.0_wp
+          rad_lw_in  = 0.0_wp
+          rad_lw_out = 0.0_wp
+
+!
+!--       Overwrite albedo if manually set in parameter file
+          IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp )  THEN
+             albedo = albedo_pars(2,albedo_type)
+          ENDIF
+    
+          alpha = albedo
+  
+!
+!--    Initialization actions for RRTMG
+       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
+#if defined ( __rrtmg )
+!
+!--       Allocate albedos
+          ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
+
+          IF ( albedo_type /= 0 )  THEN
+             IF ( albedo_lw_dif == 9999999.9_wp )  THEN
+                albedo_lw_dif = albedo_pars(0,albedo_type)
+                albedo_lw_dir = albedo_lw_dif
+             ENDIF
+             IF ( albedo_sw_dif == 9999999.9_wp )  THEN
+                albedo_sw_dif = albedo_pars(1,albedo_type)
+                albedo_sw_dir = albedo_sw_dif
+             ENDIF
+          ENDIF
+
+          aldif(:,:) = albedo_lw_dif
+          aldir(:,:) = albedo_lw_dir
+          asdif(:,:) = albedo_sw_dif
+          asdir(:,:) = albedo_sw_dir
+!
+!--       Calculate initial values of current (cosine of) the zenith angle and 
+!--       whether the sun is up
+          CALL calc_zenith     
+!
+!--       Calculate initial surface albedo
+          IF ( .NOT. constant_albedo )  THEN
+             CALL calc_albedo
+          ELSE
+             rrtm_aldif(0,:,:) = aldif(:,:)
+             rrtm_aldir(0,:,:) = aldir(:,:)
+             rrtm_asdif(0,:,:) = asdif(:,:) 
+             rrtm_asdir(0,:,:) = asdir(:,:)   
+          ENDIF
+
+!
+!--       Allocate surface emissivity
+          ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
+          rrtm_emis = emissivity
+
+!
+!--       Allocate 3d arrays of radiative fluxes and heating rates
+          ALLOCATE ( rad_sw_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( rad_sw_cs_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+          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) )
+             rad_sw_out = 0.0_wp
+          ENDIF
+
+          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
+             ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ENDIF
+
+          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
+             ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ENDIF
+
+          ALLOCATE ( rad_lw_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( rad_lw_cs_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+          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
+
+          rad_sw_hr     = 0.0_wp
+          rad_sw_cs_in  = 0.0_wp
+          rad_sw_cs_out = 0.0_wp
+          rad_sw_cs_hr  = 0.0_wp
+
+          rad_lw_hr     = 0.0_wp
+          rad_lw_cs_in  = 0.0_wp
+          rad_lw_cs_out = 0.0_wp
+          rad_lw_cs_hr  = 0.0_wp
+
+!
+!--       Allocate dummy array for storing surface temperature
+          ALLOCATE ( rrtm_tsfc(1) )
+
+!
+!--       Initialize RRTMG
+          IF ( lw_radiation )  CALL rrtmg_lw_ini ( cp )
+          IF ( sw_radiation )  CALL rrtmg_sw_ini ( cp )
+
+!
+!--       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
+
+
 !
 !--    Add timeseries for radiation model
+       dots_rad = dots_num + 1
        dots_label(dots_num+1) = "rad_net"
-       dots_label(dots_num+2) = "rad_sw_in"
-       dots_unit(dots_num+1:dots_num+2) = "W/m2"
-
-       dots_rad  = dots_num + 1
-       dots_num  = dots_num + 2
+       dots_label(dots_num+2) = "rad_lw_in"
+       dots_label(dots_num+3) = "rad_lw_out"
+       dots_label(dots_num+4) = "rad_sw_in"
+       dots_label(dots_num+5) = "rad_sw_out"
+       dots_unit(dots_num+1:dots_num+5) = "W/m2"
+       dots_num  = dots_num + 5
+
+!
+!--    Output of albedos is only required for RRTMG
+       IF ( radiation_scheme == 'rrtmg' )  THEN
+          dots_label(dots_num+1) = "rrtm_aldif"
+          dots_label(dots_num+2) = "rrtm_aldir"
+          dots_label(dots_num+3) = "rrtm_asdif"
+          dots_label(dots_num+4) = "rrtm_asdir"
+          dots_unit(dots_num+1:dots_num+4) = ""
+          dots_num  = dots_num + 4
+       ENDIF
+
+!
+!--    Calculate radiative fluxes at model start
+       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+          IF ( radiation_scheme == 'clear-sky' )  THEN
+              CALL radiation_clearsky
+          ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
+             CALL radiation_rrtmg
+          ENDIF
+       ENDIF
 
        RETURN
@@ -196,13 +599,267 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE radiation_clearsky
-    
+
+       USE indices,                                                            &
+           ONLY:  nbgp
 
        IMPLICIT NONE
 
+       INTEGER(iwp) :: i, j, k
+
+!
+!--    Calculate current zenith angle
+       CALL calc_zenith
+
+!
+!--    Calculate sky transmissivity
+       sky_trans = 0.6_wp + 0.2_wp * zenith(0)
+
+!
+!--    Calculate value of the Exner function
+       exn = (surface_pressure / 1000.0_wp )**0.286_wp
+!
+!--    Calculate radiation fluxes and net radiation (rad_net) for each grid 
+!--    point
+       DO i = nxl, nxr
+          DO j = nys, nyn
+             k = nzb_s_inner(j,i)
+             rad_sw_in(0,j,i)  = solar_constant * sky_trans * zenith(0)
+             rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
+             rad_lw_out(0,j,i) = sigma_sb * (pt(k,j,i) * exn)**4
+             rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn)**4
+             rad_net(j,i)      = rad_sw_in(0,j,i) - rad_sw_out(0,j,i)          &
+                                 + rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
+
+          ENDDO
+       ENDDO
+
+       CALL exchange_horiz_2d( rad_lw_in,  nbgp )
+       CALL exchange_horiz_2d( rad_lw_out, nbgp )
+       CALL exchange_horiz_2d( rad_sw_in,  nbgp )
+       CALL exchange_horiz_2d( rad_sw_out, nbgp )
+       CALL exchange_horiz_2d( rad_net,    nbgp )
+
+       RETURN
+
+    END SUBROUTINE radiation_clearsky
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!-- Implementation of the RRTMG radiation_scheme
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_rrtmg
+
+       USE indices,                                                            &
+           ONLY:  nbgp
+
+       USE particle_attributes,                                                &
+           ONLY:  grid_particles, number_of_particles, particles,              &
+                  particle_advection_start, prt_count
+
+       IMPLICIT NONE
+
+#if defined ( __rrtmg )
+
+       INTEGER(iwp) :: i, j, k, n          !: 
+
+       REAL(wp)     ::  s_r2      !: weighted sum over all droplets with r^2
+       REAL(wp)     ::  s_r3      !: weighted sum over all droplets with r^3
+
+!
+!--    Calculate current (cosine of) zenith angle and whether the sun is up
+       CALL calc_zenith     
+!
+!--    Calculate surface albedo
+       IF ( .NOT. constant_albedo )  THEN
+          CALL calc_albedo
+       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
+!
+!--    Loop over all grid points
+       DO i = nxl, nxr
+          DO j = nys, nyn
+
+!
+!--          Prepare profiles of temperature and H2O volume mixing ratio
+             rrtm_tlev(0,nzb+1) = pt(nzb,j,i)
+
+             DO k = nzb+1, nzt+1
+                rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp         &
+                                 )**0.286_wp + ( l_v / cp ) * ql(k,j,i)
+                rrtm_h2ovmr(0,k) = mol_weight_air_d_wv * (q(k,j,i) - ql(k,j,i))
+
+             ENDDO
+
+!
+!--          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/mÂ² instead of kg/kg m.
+             rrtm_cldfr  = 0.0_wp
+             rrtm_reliq  = 0.0_wp
+             rrtm_cliqwp = 0.0_wp
+
+             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) .GT. 0 )  THEN
+                   rrtm_cldfr(0,k) = 1.0_wp
+                   rrtm_icld = 1
+
+!
+!--                Calculate cloud droplet effective radius
+                   IF ( cloud_physics )  THEN
+                      rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)          &
+                                      / ( 4.0_wp * pi * nc_const * 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) .GT. 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
+                ELSE
+                   rrtm_icld = 0
+                ENDIF
+             ENDDO
+
+!
+!--          Set surface temperature
+             rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
+
+             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 )
+
+                DO k = nzb, nzt+1
+                   rad_lw_in(k,j,i)  = rrtm_lwdflx(0,k)
+                   rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
+                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(:,j,i),&
+               rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
+               0.0_wp          , day          , 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 )
+  
+                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
+             ENDIF
+
+!
+!--          Calculate surface net radiation
+             rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i)           &
+                            + rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
+
+          ENDDO
+       ENDDO
+
+       CALL exchange_horiz( rad_lw_in,  nbgp )
+       CALL exchange_horiz( rad_lw_out, nbgp )
+       CALL exchange_horiz( rad_sw_in,  nbgp )
+       CALL exchange_horiz( rad_sw_out, nbgp ) 
+       CALL exchange_horiz_2d( rad_net, nbgp )
+#endif
+
+    END SUBROUTINE radiation_rrtmg
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!-- Calculate the cosine of the zenith angle (variable is called zenith)
+!------------------------------------------------------------------------------!
+    SUBROUTINE calc_zenith
+
+       IMPLICIT NONE
+
+       REAL(wp) ::  declination,  & !: solar declination angle
+                    hour_angle      !: solar hour angle
 !
 !--    Calculate current day and time based on the initial values and simulation
 !--    time
-       day = day_init + FLOOR( (time_utc_init + time_since_reference_point)    &
-                               / 86400.0_wp )
+       day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point)    &
+                               / 86400.0_wp ), KIND=iwp)
        time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
 
@@ -210,62 +867,720 @@
 !
 !--    Calculate solar declination and hour angle   
-       declination = ASIN( decl_1 * SIN(decl_2 * day - decl_3) )
+       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
        hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
 
 !
 !--    Calculate zenith angle
-       zenith = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)       &
+       zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)      &
                                             * COS(hour_angle)
-       zenith = MAX(0.0_wp,zenith)
-
-!
-!--    Calculate sky transmissivity
-       sky_trans = 0.6_wp + 0.2_wp * zenith
-
-!
-!--    Calculate value of the Exner function
-       exn = (surface_pressure / 1000.0_wp )**0.286_wp
-
-!
-!--    Calculate radiation fluxes and net radiation (rad_net) for each grid 
-!--    point
-       DO i = nxlg, nxrg
-          DO j = nysg, nyng
-
-             k = nzb_s_inner(j,i)
-             rad_sw_in(j,i)  = sw_0 * sky_trans * zenith
-             rad_sw_out(j,i) = - alpha(j,i) * rad_sw_in(j,i)
-             rad_lw_out(j,i) = - sigma_sb * (pt(k,j,i) * exn)**4
-             rad_lw_in(j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn)**4
-             rad_net(j,i)    = rad_sw_in(j,i) + rad_sw_out(j,i)                &
-                                + rad_lw_in(j,i) + rad_lw_out(j,i)
-
+       zenith(0) = MAX(0.0_wp,zenith(0))
+
+!
+!--    Check if the sun is up (otheriwse shortwave calculations can be skipped)
+       IF ( zenith(0) .GT. 0.0_wp )  THEN
+          sun_up = .TRUE.
+       ELSE
+          sun_up = .FALSE.
+       END IF
+
+    END SUBROUTINE calc_zenith
+
+#if defined ( __rrtmg )
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!-- Calculates surface albedo components based on Briegleb (1992) and 
+!-- Briegleb et al. (1986)
+!------------------------------------------------------------------------------!
+    SUBROUTINE calc_albedo
+
+        IMPLICIT NONE
+
+        IF ( sun_up )  THEN
+!
+!--        Ocean
+           IF ( albedo_type == 1 )  THEN
+              rrtm_aldir(0,:,:) = 0.026_wp / ( zenith(0)**1.7_wp + 0.065_wp )  &
+                                  + 0.15_wp * ( zenith(0) - 0.1_wp )           &
+                                            * ( zenith(0) - 0.5_wp )           &
+                                            * ( zenith(0) - 1.0_wp )
+              rrtm_asdir(0,:,:) = rrtm_aldir(0,:,:)
+!
+!--        Snow
+           ELSEIF ( albedo_type == 16 )  THEN
+              IF ( zenith(0) < 0.5_wp )  THEN
+                 rrtm_aldir(0,:,:) = 0.5_wp * (1.0_wp - aldif)                 &
+                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
+                                     * zenith(0))) - 1.0_wp
+                 rrtm_asdir(0,:,:) = 0.5_wp * (1.0_wp - asdif)                 &
+                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
+                                     * zenith(0))) - 1.0_wp
+
+                 rrtm_aldir(0,:,:) = MIN(0.98_wp, rrtm_aldir(0,:,:))
+                 rrtm_asdir(0,:,:) = MIN(0.98_wp, rrtm_asdir(0,:,:))
+              ELSE
+                 rrtm_aldir(0,:,:) = aldif
+                 rrtm_asdir(0,:,:) = asdif
+              ENDIF
+!
+!--        Sea ice
+           ELSEIF ( albedo_type == 15 )  THEN
+                 rrtm_aldir(0,:,:) = aldif
+                 rrtm_asdir(0,:,:) = asdif
+!
+!--        Land surfaces
+           ELSE
+              SELECT CASE ( albedo_type )
+
+!
+!--              Surface types with strong zenith dependence
+                 CASE ( 1, 2, 3, 4, 11, 12, 13 )
+                    rrtm_aldir(0,:,:) = aldif * 1.4_wp /                       &
+                                        (1.0_wp + 0.8_wp * zenith(0))
+                    rrtm_asdir(0,:,:) = asdif * 1.4_wp /                       &
+                                        (1.0_wp + 0.8_wp * zenith(0))
+!
+!--              Surface types with weak zenith dependence
+                 CASE ( 5, 6, 7, 8, 9, 10, 14 )
+                    rrtm_aldir(0,:,:) = aldif * 1.1_wp /                       &
+                                        (1.0_wp + 0.2_wp * zenith(0))
+                    rrtm_asdir(0,:,:) = asdif * 1.1_wp /                       &
+                                        (1.0_wp + 0.2_wp * zenith(0))
+
+                 CASE DEFAULT
+
+              END SELECT
+           ENDIF
+!
+!--        Diffusive albedo is taken from Table 2
+           rrtm_aldif(0,:,:) = aldif
+           rrtm_asdif(0,:,:) = asdif
+
+        ELSE
+
+           rrtm_aldir(0,:,:) = 0.0_wp
+           rrtm_asdir(0,:,:) = 0.0_wp
+           rrtm_aldif(0,:,:) = 0.0_wp
+           rrtm_asdif(0,:,:) = 0.0_wp
+        ENDIF
+    END SUBROUTINE calc_albedo
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!-- Read sounding data (pressure and temperature) from RADIATION_DATA.
+!------------------------------------------------------------------------------!
+    SUBROUTINE read_sounding_data
+
+       USE netcdf_control
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) :: id, id_dim_zrad, id_var, k, nz_snd, nz_snd_start, nz_snd_end
+
+       REAL(wp) :: t_surface
+
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  hyp_snd_tmp, t_snd_tmp
+
+!
+!--    In case of updates, deallocate arrays first (sufficient to check one
+!--    array as the others are automatically allocated). This is required
+!--    because nzt_rad might change during the update
+       IF ( ALLOCATED ( hyp_snd ) )  THEN
+          DEALLOCATE( hyp_snd )
+          DEALLOCATE( t_snd )
+          DEALLOCATE( q_snd  )
+          DEALLOCATE ( rrtm_play )
+          DEALLOCATE ( rrtm_plev )
+          DEALLOCATE ( rrtm_tlay )
+          DEALLOCATE ( rrtm_tlev )
+          DEALLOCATE ( rrtm_h2ovmr )
+          DEALLOCATE ( rrtm_cicewp )
+          DEALLOCATE ( rrtm_cldfr )
+          DEALLOCATE ( rrtm_cliqwp )
+          DEALLOCATE ( rrtm_reice )
+          DEALLOCATE ( rrtm_reliq )
+          DEALLOCATE ( rrtm_lw_taucld )
+          DEALLOCATE ( rrtm_lw_tauaer )
+          DEALLOCATE ( rrtm_lwdflx  )
+          DEALLOCATE ( rrtm_lwuflx  )
+          DEALLOCATE ( rrtm_lwhr  )
+          DEALLOCATE ( rrtm_lwuflxc )
+          DEALLOCATE ( rrtm_lwdflxc )
+          DEALLOCATE ( rrtm_lwhrc )
+          DEALLOCATE ( rrtm_sw_taucld )
+          DEALLOCATE ( rrtm_sw_ssacld )
+          DEALLOCATE ( rrtm_sw_asmcld )
+          DEALLOCATE ( rrtm_sw_fsfcld )
+          DEALLOCATE ( rrtm_sw_tauaer )
+          DEALLOCATE ( rrtm_sw_ssaaer )
+          DEALLOCATE ( rrtm_sw_asmaer ) 
+          DEALLOCATE ( rrtm_sw_ecaer )    
+          DEALLOCATE ( rrtm_swdflx  )
+          DEALLOCATE ( rrtm_swuflx  )
+          DEALLOCATE ( rrtm_swhr  )
+          DEALLOCATE ( rrtm_swuflxc )
+          DEALLOCATE ( rrtm_swdflxc )
+          DEALLOCATE ( rrtm_swhrc )
+       ENDIF
+
+!
+!--    Open file for reading
+       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
+       CALL handle_netcdf_error( 'netcdf', 549 )
+
+!
+!--    Inquire dimension of z axis and save in nz_snd
+       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
+       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
+       CALL handle_netcdf_error( 'netcdf', 551 )
+
+!
+! !--    Allocate temporary array for storing pressure data
+       ALLOCATE( hyp_snd_tmp(nzb+1:nz_snd) )
+       hyp_snd_tmp = 0.0_wp
+
+
+!--    Read pressure from file
+       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
+       nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/),    &
+                               count = (/nz_snd/) )
+       CALL handle_netcdf_error( 'netcdf', 552 )
+
+!
+!--    Allocate temporary array for storing temperature data
+       ALLOCATE( t_snd_tmp(nzb+1:nz_snd) )
+       t_snd_tmp = 0.0_wp
+
+!
+!--    Read temperature from file
+       nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
+       nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/),      &
+                               count = (/nz_snd/) )
+       CALL handle_netcdf_error( 'netcdf', 553 )
+
+!
+!--    Calculate start of sounding data
+       nz_snd_start = nz_snd + 1
+       nz_snd_end   = nz_snd_end
+
+!
+!--    Start filling vertical dimension at 10hPa above the model domain (hyp is
+!--    in Pa, hyp_snd in hPa).
+       DO  k = 1, nz_snd
+          IF ( hyp_snd_tmp(k) .LT. (hyp(nzt+1) - 1000.0_wp) * 0.01_wp )  THEN
+             nz_snd_start = k
+             EXIT
+          END IF
+       END DO
+
+       IF ( nz_snd_start .LE. nz_snd )  THEN
+          nz_snd_end = nz_snd - 1
+       END IF
+
+
+!
+!--    Calculate of total grid points for RRTMG calculations
+       nzt_rad = nzt + nz_snd_end - nz_snd_start + 2
+
+!
+!--    Save data above LES domain in hyp_snd, t_snd and q_snd
+!--    Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
+       ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
+       ALLOCATE( t_snd(nzb+1:nzt_rad)   )
+       ALLOCATE( q_snd(nzb+1:nzt_rad)   )
+       hyp_snd = 0.0_wp
+       t_snd = 0.0_wp
+       q_snd = 0.0_wp
+
+       hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start:nz_snd_end)
+       t_snd(nzt+2:nzt_rad)   = t_snd_tmp(nz_snd_start:nz_snd_end)
+
+       DEALLOCATE ( hyp_snd_tmp )
+       DEALLOCATE ( t_snd_tmp )
+
+       nc_stat = NF90_CLOSE( id )
+
+!
+!--    Calculate pressure levels on zu and zw grid. Sounding data is added at 
+!--    top of the LES domain. This routine does not consider horizontal or 
+!--    vertical variability of pressure and temperature
+       ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1)   )
+       ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2)   )
+
+       t_surface = pt_surface * (surface_pressure / 1000.0_wp )**0.286_wp
+       DO k = nzb+1, nzt+1
+          rrtm_play(0,k) = hyp(k) * 0.01_wp
+          rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
+                         t_surface )**(1.0_wp/0.286_wp)
+       ENDDO
+
+       DO k = nzt+2, nzt_rad
+          rrtm_play(0,k) = hyp_snd(k)
+          rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
+       ENDDO
+       rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad),                   &
+                                   1.5 * hyp_snd(nzt_rad)                      &
+                                 - 0.5 * hyp_snd(nzt_rad-1) )
+       rrtm_plev(0,nzt_rad+2)  = MIN( 1.0E-4_wp,                               &
+                                      0.25_wp * rrtm_plev(0,nzt_rad+1) )
+
+       rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
+
+!
+!--    Calculate temperature/humidity levels at top of the LES domain. 
+!--    Currently, the temperature is taken from sounding data (might lead to a 
+!--    temperature jump at interface. To do: Humidity is currently not 
+!--    calculated above the LES domain.
+       ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1)   )
+       ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2)   )
+       ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
+
+       DO k = nzt+8, nzt_rad
+          rrtm_tlay(0,k)   = t_snd(k)
+          rrtm_h2ovmr(0,k) = q_snd(k)
+       ENDDO
+       rrtm_tlay(0,nzt_rad+1)   = 2.0_wp * rrtm_tlay(0,nzt_rad)                &
+                                  - rrtm_tlay(0,nzt_rad-1)
+       DO k = nzt+9, nzt_rad+1
+          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
+       rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
+
+       rrtm_tlev(0,nzt_rad+2)   = 2.0_wp * rrtm_tlay(0,nzt_rad+1)              &
+                                  - rrtm_tlev(0,nzt_rad)
+!
+!--    Allocate remaining RRTMG arrays
+       ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
+       ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
+       ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
+       ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
+       ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
+       ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
+       ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
+       ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
+       ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
+       ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
+       ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
+       ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
+       ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
+       ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) ) 
+       ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )    
+
+!
+!--    The ice phase is currently not considered in PALM
+       rrtm_cicewp = 0.0_wp
+       rrtm_reice  = 0.0_wp
+
+!
+!--    Set other parameters (move to NAMELIST parameters in the future)
+       rrtm_lw_tauaer = 0.0_wp
+       rrtm_lw_taucld = 0.0_wp
+       rrtm_sw_taucld = 0.0_wp
+       rrtm_sw_ssacld = 0.0_wp
+       rrtm_sw_asmcld = 0.0_wp
+       rrtm_sw_fsfcld = 0.0_wp
+       rrtm_sw_tauaer = 0.0_wp
+       rrtm_sw_ssaaer = 0.0_wp
+       rrtm_sw_asmaer = 0.0_wp
+       rrtm_sw_ecaer  = 0.0_wp
+
+
+       ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1)  )
+       ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1)  )
+       ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1)  )
+       ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
+       ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
+       ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
+
+       rrtm_swdflx  = 0.0_wp
+       rrtm_swuflx  = 0.0_wp
+       rrtm_swhr    = 0.0_wp  
+       rrtm_swuflxc = 0.0_wp
+       rrtm_swdflxc = 0.0_wp
+       rrtm_swhrc   = 0.0_wp
+
+       ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1)  )
+       ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1)  )
+       ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1)  )
+       ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
+       ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
+       ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
+
+       rrtm_lwdflx  = 0.0_wp
+       rrtm_lwuflx  = 0.0_wp
+       rrtm_lwhr    = 0.0_wp  
+       rrtm_lwuflxc = 0.0_wp
+       rrtm_lwdflxc = 0.0_wp
+       rrtm_lwhrc   = 0.0_wp
+
+
+    END SUBROUTINE read_sounding_data
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!-- Read trace gas data from file
+!------------------------------------------------------------------------------!
+    SUBROUTINE read_trace_gas_data
+
+       USE netcdf_control
+       USE rrsw_ncpar
+
+       IMPLICIT NONE
+
+       INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !: number of trace gases
+
+       CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER ::              &
+           trace_names = (/'O3   ', 'CO2  ', 'CH4  ', 'N2O  ', 'O2   ',        &
+                           'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
+
+       INTEGER(iwp) :: id, k, m, n, nabs, np, id_abs, id_dim, id_var
+
+       REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
+
+
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  p_mls,         & !: pressure levels for the absorbers
+                                                 rrtm_play_tmp, & !: temporary array for pressure zu-levels
+                                                 rrtm_plev_tmp, & !: temporary array for pressure zw-levels
+                                                 trace_path_tmp   !: temporary array for storing trace gas path data
+
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  trace_mls,      & !: array for storing the absorber amounts
+                                                 trace_mls_path, & !: array for storing trace gas path data
+                                                 trace_mls_tmp     !: temporary array for storing trace gas data
+
+
+!
+!--    In case of updates, deallocate arrays first (sufficient to check one
+!--    array as the others are automatically allocated)
+       IF ( ALLOCATED ( rrtm_o3vmr ) )  THEN
+          DEALLOCATE ( rrtm_o3vmr  )
+          DEALLOCATE ( rrtm_co2vmr )
+          DEALLOCATE ( rrtm_ch4vmr )
+          DEALLOCATE ( rrtm_n2ovmr )
+          DEALLOCATE ( rrtm_o2vmr  )
+          DEALLOCATE ( rrtm_cfc11vmr )
+          DEALLOCATE ( rrtm_cfc12vmr )
+          DEALLOCATE ( rrtm_cfc22vmr )
+          DEALLOCATE ( rrtm_ccl4vmr  )
+       ENDIF
+
+!
+!--    Allocate trace gas profiles
+       ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1)  )
+       ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
+       ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
+       ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
+       ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1)  )
+       ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
+       ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
+       ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
+       ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1)  )
+
+!
+!--    Open file for reading
+       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
+       CALL handle_netcdf_error( 'netcdf', 549 )
+!
+!--    Inquire dimension ids and dimensions
+       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
+       CALL handle_netcdf_error( 'netcdf', 550 )
+       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np) 
+       CALL handle_netcdf_error( 'netcdf', 550 )
+
+       nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
+       CALL handle_netcdf_error( 'netcdf', 550 )
+       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs ) 
+       CALL handle_netcdf_error( 'netcdf', 550 )
+   
+
+!
+!--    Allocate pressure, and trace gas arrays     
+       ALLOCATE( p_mls(1:np) )
+       ALLOCATE( trace_mls(1:num_trace_gases,1:np) ) 
+       ALLOCATE( trace_mls_tmp(1:nabs,1:np) ) 
+
+
+       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
+       CALL handle_netcdf_error( 'netcdf', 550 )
+       nc_stat = NF90_GET_VAR( id, id_var, p_mls )
+       CALL handle_netcdf_error( 'netcdf', 550 )
+
+       nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
+       CALL handle_netcdf_error( 'netcdf', 550 )
+       nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
+       CALL handle_netcdf_error( 'netcdf', 550 )
+
+
+!
+!--    Write absorber amounts (mls) to trace_mls
+       DO n = 1, num_trace_gases
+          CALL getAbsorberIndex( TRIM( trace_names(n) ), id_abs )
+
+          trace_mls(n,1:np) = trace_mls_tmp(id_abs,1:np)
+
+!
+!--       Replace missing values by zero
+          WHERE ( trace_mls(n,:) > 2.0_wp )  
+             trace_mls(n,:) = 0.0_wp
+          END WHERE
+       END DO
+
+       DEALLOCATE ( trace_mls_tmp )
+
+       nc_stat = NF90_CLOSE( id )
+       CALL handle_netcdf_error( 'netcdf', 551 )
+
+!
+!--    Add extra pressure level for calculations of the trace gas paths
+       ALLOCATE ( rrtm_play_tmp(1:nzt_rad+1) )
+       ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
+
+       rrtm_play_tmp(1:nzt_rad)   = rrtm_play(0,1:nzt_rad) 
+       rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
+       rrtm_play_tmp(nzt_rad+1)   = rrtm_plev(0,nzt_rad+1) * 0.5_wp
+       rrtm_plev_tmp(nzt_rad+2)   = MIN( 1.0E-4_wp, 0.25_wp                    &
+                                         * rrtm_plev(0,nzt_rad+1) )
+ 
+!
+!--    Calculate trace gas path (zero at surface) with interpolation to the
+!--    sounding levels
+       ALLOCATE ( trace_mls_path(1:nzt_rad+2,1:num_trace_gases) )
+
+       trace_mls_path(nzb+1,:) = 0.0_wp
+       
+       DO k = nzb+2, nzt_rad+2
+          DO m = 1, num_trace_gases
+             trace_mls_path(k,m) = trace_mls_path(k-1,m)
+
+!
+!--          When the pressure level is higher than the trace gas pressure
+!--          level, assume that 
+             IF ( rrtm_plev_tmp(k-1) .GT. p_mls(1) )  THEN             
+                
+                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1)     &
+                                      * ( rrtm_plev_tmp(k-1)                   &
+                                          - MAX( p_mls(1), rrtm_plev_tmp(k) )  &
+                                        ) / g
+             ENDIF
+
+!
+!--          Integrate for each sounding level from the contributing p_mls 
+!--          levels
+             DO n = 2, np
+!
+!--             Limit p_mls so that it is within the model level
+                p_mls_u = MIN( rrtm_plev_tmp(k-1),                             &
+                          MAX( rrtm_plev_tmp(k), p_mls(n) ) )
+                p_mls_l = MIN( rrtm_plev_tmp(k-1),                             &
+                          MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
+
+                IF ( p_mls_l .GT. p_mls_u )  THEN
+
+!
+!--                Calculate weights for interpolation
+                   p_mls_m = 0.5_wp * (p_mls_l + p_mls_u)
+                   p_wgt_u = (p_mls(n-1) - p_mls_m) / (p_mls(n-1) - p_mls(n))
+                   p_wgt_l = (p_mls_m - p_mls(n))   / (p_mls(n-1) - p_mls(n))
+
+!
+!--                Add level to trace gas path
+                   trace_mls_path(k,m) = trace_mls_path(k,m)                   &
+                                         +  ( p_wgt_u * trace_mls(m,n)         &
+                                            + p_wgt_l * trace_mls(m,n-1) )     &
+                                         * (p_mls_l + p_mls_u) / g
+                ENDIF
+             ENDDO
+
+             IF ( rrtm_plev_tmp(k) .LT. p_mls(np) )  THEN 
+                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np)    &
+                                      * ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
+                                          - rrtm_plev_tmp(k)                   &
+                                        ) / g 
+             ENDIF  
           ENDDO
        ENDDO
 
-       RETURN
-
-    END SUBROUTINE radiation_clearsky
+
+!
+!--    Prepare trace gas path profiles
+       ALLOCATE ( trace_path_tmp(1:nzt_rad+1) )
+
+       DO m = 1, num_trace_gases
+
+          trace_path_tmp(1:nzt_rad+1) = ( trace_mls_path(2:nzt_rad+2,m)        &
+                                       - trace_mls_path(1:nzt_rad+1,m) ) * g   &
+                                       / ( rrtm_plev_tmp(1:nzt_rad+1)          &
+                                       - rrtm_plev_tmp(2:nzt_rad+2) )
+
+!
+!--       Save trace gas paths to the respective arrays
+          SELECT CASE ( TRIM( trace_names(m) ) )
+
+             CASE ( 'O3' )
+
+                rrtm_o3vmr(0,:) = trace_path_tmp(:)
+
+             CASE ( 'CO2' )
+
+                rrtm_co2vmr(0,:) = trace_path_tmp(:)
+
+             CASE ( 'CH4' )
+
+                rrtm_ch4vmr(0,:) = trace_path_tmp(:)
+
+             CASE ( 'N2O' )
+
+                rrtm_n2ovmr(0,:) = trace_path_tmp(:)
+
+             CASE ( 'O2' )
+
+                rrtm_o2vmr(0,:) = trace_path_tmp(:)
+
+             CASE ( 'CFC11' )
+
+                rrtm_cfc11vmr(0,:) = trace_path_tmp(:)
+
+             CASE ( 'CFC12' )
+
+                rrtm_cfc12vmr(0,:) = trace_path_tmp(:)
+
+             CASE ( 'CFC22' )
+
+                rrtm_cfc22vmr(0,:) = trace_path_tmp(:)
+
+             CASE ( 'CCL4' )
+
+                rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
+
+             CASE DEFAULT
+
+          END SELECT
+
+       ENDDO
+
+       DEALLOCATE ( trace_path_tmp )
+       DEALLOCATE ( trace_mls_path )
+       DEALLOCATE ( rrtm_play_tmp )
+       DEALLOCATE ( rrtm_plev_tmp )
+       DEALLOCATE ( trace_mls )
+       DEALLOCATE ( p_mls )
+
+    END SUBROUTINE read_trace_gas_data
+
+#endif
+
 
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!-- Prescribed net radiation
-!------------------------------------------------------------------------------!
-    SUBROUTINE radiation_constant
-
-       rad_net = net_radiation
-
-    END SUBROUTINE radiation_constant
+!-- Calculate temperature tendency due to radiative cooling/heating.
+!-- Cache-optimized version.
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_tendency_ij ( i, j, tend )
+
+       USE arrays_3d,                                                          &
+           ONLY:  dzw
+
+       USE cloud_parameters,                                                   &
+           ONLY:  pt_d_t, cp
+
+       USE control_parameters,                                                 &
+           ONLY:  rho_surface
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) :: i, j, k
+
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend
+
+#if defined ( __rrtmg )
+
+       REAL(wp) :: rad_sw_net_l, rad_sw_net_u, rad_lw_net_l, rad_lw_net_u
+
+       rad_sw_net_l = 0.0_wp
+       rad_sw_net_u = 0.0_wp
+       rad_lw_net_l = 0.0_wp
+       rad_lw_net_u = 0.0_wp
+
+!
+!--    Calculate radiative flux divergence
+       DO k = nzb+1, nzt+1
+
+          rad_sw_net_l = rad_sw_in(k-1,j,i) - rad_sw_out(k-1,j,i)
+          rad_sw_net_u = rad_sw_in(k,j,i)   - rad_sw_out(k,j,i)
+          rad_lw_net_l = rad_lw_in(k-1,j,i) - rad_lw_out(k-1,j,i)
+          rad_lw_net_u = rad_lw_in(k,j,i)   - rad_lw_out(k,j,i)
+
+          tend(k,j,i) = tend(k,j,i) + ( rad_sw_net_u - rad_sw_net_l            &
+                                      + rad_lw_net_u - rad_lw_net_l ) /        &
+                                     ( dzw(k) * cp * rho_surface ) * pt_d_t(k) 
+       ENDDO
+
+#endif
+
+    END SUBROUTINE radiation_tendency_ij
+
 
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!-- Implementation of the RRTM radiation_scheme
-!------------------------------------------------------------------------------!
-    SUBROUTINE radiation_rrtm
-
-
-    END SUBROUTINE radiation_rrtm
+!-- Calculate temperature tendency due to radiative cooling/heating.
+!-- Vector-optimized version
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_tendency ( tend )
+
+       USE arrays_3d,                                                          &
+           ONLY:  dzw
+
+       USE cloud_parameters,                                                   &
+           ONLY:  pt_d_t, cp
+
+       USE indices,                                                            &
+           ONLY:  nxl, nxr, nyn, nys
+
+       USE control_parameters,                                                 &
+           ONLY:  rho_surface
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) :: i, j, k
+
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend
+
+#if defined ( __rrtmg )
+
+       REAL(wp) :: rad_sw_net_l, rad_sw_net_u, rad_lw_net_l, rad_lw_net_u
+
+       rad_sw_net_l = 0.0_wp
+       rad_sw_net_u = 0.0_wp
+       rad_lw_net_l = 0.0_wp
+       rad_lw_net_u = 0.0_wp
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+
+!
+!--          Calculate radiative flux divergence
+             DO k = nzb+1, nzt+1
+
+                rad_sw_net_l = rad_sw_in(k-1,j,i) - rad_sw_out(k-1,j,i)
+                rad_sw_net_u = rad_sw_in(k  ,j,i) - rad_sw_out(k  ,j,i)
+                rad_lw_net_l = rad_lw_in(k-1,j,i) - rad_lw_out(k-1,j,i)
+                rad_lw_net_u = rad_lw_in(k  ,j,i) - rad_lw_out(k  ,j,i)
+
+                tend(k,j,i) = tend(k,j,i) + ( rad_sw_net_u - rad_sw_net_l      &
+                                      + rad_lw_net_u - rad_lw_net_l ) /        &
+                                      ( dzw(k) * cp * rho_surface ) * pt_d_t(k) 
+             ENDDO
+         ENDDO
+       ENDDO
+#endif
+
+    END SUBROUTINE radiation_tendency
 
  END MODULE radiation_model_mod
