Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 1584)
+++ palm/trunk/SOURCE/Makefile	(revision 1585)
@@ -20,5 +20,5 @@
 # Current revisions:
 # ------------------
-# 
+# Added user_init_radiation.f90
 # 
 # Former revisions:
@@ -223,5 +223,5 @@
 	init_pegrid.f90 init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \
 	interaction_droplets_ptq.f90 land_surface_model.f90 local_flush.f90 \
-   local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \
+	local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \
 	local_tremain_ini.f90 lpm.f90 lpm_advec.f90 lpm_boundary_conds.f90 \
 	lpm_calc_liquid_water_content.f90 lpm_collision_kernels.f90 \
@@ -253,9 +253,11 @@
 	user_dvrp_coltab.f90 user_header.f90 user_init.f90 \
 	user_init_3d_model.f90 user_init_grid.f90 user_init_land_surface.f90 \
-	user_init_plant_canopy.f90 user_last_actions.f90 user_lpm_advec.f90 \
+	user_init_plant_canopy.f90 user_init_radiation.f90 \
+	user_last_actions.f90 user_lpm_advec.f90 \
 	user_lpm_init.f90 user_lpm_set_attributes.f90 user_module.f90 \
 	user_parin.f90 user_read_restart_data.f90 \
 	user_spectra.f90 user_statistics.f90 wall_fluxes.f90 \
 	write_3d_binary.f90 write_var_list.f90
+
 
 OBJS=$(SOURCES:.f90=.o)
@@ -427,5 +429,5 @@
 	nudging.o plant_canopy_model.o production_e.o subsidence.o user_actions.o
 progress_bar.o: modules.o mod_kinds.o
-radiation_model.o : modules.o
+radiation_model.o : modules.o mod_particle_attributes.o
 random_function.o: mod_kinds.o
 random_gauss.o: mod_kinds.o random_function.o random_generator_parallel.o
@@ -473,4 +475,5 @@
 user_init_land_surface.o: modules.o mod_kinds.o user_module.o land_surface_model.o
 user_init_plant_canopy.o: modules.o mod_kinds.o user_module.o plant_canopy_model.o
+user_init_radiation.o: modules.o mod_kinds.o user_module.o radiation_model.o
 user_last_actions.o: modules.o mod_kinds.o user_module.o
 user_lpm_advec.o: modules.o mod_kinds.o user_module.o
Index: palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/average_3d_data.f90	(revision 1584)
+++ palm/trunk/SOURCE/average_3d_data.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adapted for RRTMG
 ! 
 ! Former revisions:
@@ -88,5 +88,7 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av
+        ONLY:  rad_net, rad_net_av, 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
 
     IMPLICIT NONE
@@ -349,15 +351,44 @@
              ENDDO
 
-         CASE ( 'rad_sw_in*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   rad_sw_in_av(j,i) = rad_sw_in_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
          CASE ( 'rad_net*' )
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
                    rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_lw_in' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_lw_out' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_in' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_out' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
                 ENDDO
              ENDDO
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 1584)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added support for RRTMG
 !
 ! Former revisions:
@@ -1143,14 +1143,24 @@
 
     IF ( radiation )  THEN
-       IF ( radiation_scheme == 'constant' )  THEN
-          irad_scheme = 0
-       ELSEIF ( radiation_scheme == 'clear-sky' )  THEN
-          irad_scheme = 1
-       ELSEIF ( radiation_scheme == 'rrtm' )  THEN
-          irad_scheme = 2
-       ELSE
+       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
+       ENDIF
+       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
     ENDIF
@@ -3058,4 +3068,55 @@
              ENDIF
 
+          CASE ( 'rad_lw_in' )
+             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
+                                 'lable for radiation = .FALSE. or ' //        &
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 102
+                dopr_unit(i)  = 'W/m2'
+                hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'rad_lw_out' )
+             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
+                                 'lable for radiation = .FALSE. or ' //        &
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 103
+                dopr_unit(i)  = 'W/m2'
+                hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'rad_sw_in' )
+             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
+                                 'lable for radiation = .FALSE. or ' //        &
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 104
+                dopr_unit(i)  = 'W/m2'
+                hom(:,2,104,:)  = SPREAD( zw, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'rad_sw_out')
+             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
+                                 'lable for radiation = .FALSE. or ' //        &
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 105
+                dopr_unit(i)  = 'W/m2'
+                hom(:,2,105,:)  = SPREAD( zw, 2, statistic_regions+1 )
+             ENDIF
 
           CASE DEFAULT
@@ -3257,4 +3318,14 @@
              unit = 'kg/kg'
 
+
+          CASE ( 'rad_lw_in', 'rad_lw_out', '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 = 'W/m2'
+
           CASE ( 'rho' )
              IF ( .NOT. ocean )  THEN
@@ -3293,6 +3364,6 @@
                  'm_liq_eb*', 'pra*', 'prr*', 'qsws*', 'qsws_eb*',             &
                  'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*', 'rad_net*',  &
-                 'rad_sw_in*', 'r_a*', 'r_s*', 'shf*', 'shf_eb*', 't*',        &
-                 'u*', 'z0*', 'z0h*' )
+                 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*',   &
+                 'r_a*', 'r_s*', 'shf*', 'shf_eb*', 't*', 'u*', 'z0*', 'z0h*' )
              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
                 message_string = 'illegal value for data_output: "' //         &
@@ -3301,4 +3372,17 @@
                 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 ) == 'c_liq*'  .AND.  .NOT. land_surface )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -3405,6 +3489,9 @@
              IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
              IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
-             IF ( TRIM( var ) == 'rad_net*')       unit = 'W/m2'     
-             IF ( TRIM( var ) == 'rad_sw_in*')     unit = 'W/m2' 
+             IF ( TRIM( var ) == 'rad_net*'      ) 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 = '' 
              IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
              IF ( TRIM( var ) == 'r_s*')     unit = 's/m' 
Index: palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 1584)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added support for RRTMG
 ! 
 ! Former revisions:
@@ -154,5 +154,7 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av
+        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av, rad_sw_out,       &
+               rad_sw_out_av, rad_lw_in, rad_lw_in_av, rad_lw_out,             &
+               rad_lw_out_av 
 
     IMPLICIT NONE
@@ -913,21 +915,32 @@
                 level_z(nzb+1) = zu(nzb+1)
 
-             CASE ( 'rad_sw_in*_xy' )        ! 2d-array
-                IF ( av == 0 ) THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) =  rad_sw_in(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng 
-                         local_pf(i,j,nzb+1) =  rad_sw_in_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => rad_lw_in
+                ELSE
+                   to_be_resorted => rad_lw_in_av
+                ENDIF
+
+             CASE ( 'rad_lw_out_xy', 'rad_lw_out_xz', 'rad_lw_out_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => rad_lw_out
+                ELSE
+                   to_be_resorted => rad_lw_out_av
+                ENDIF
+
+             CASE ( 'rad_sw_in_xy', 'rad_sw_in_xz', 'rad_sw_in_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => rad_sw_in
+                ELSE
+                   to_be_resorted => rad_sw_in_av
+                ENDIF
+
+             CASE ( 'rad_sw_out_xy', 'rad_sw_out_xz', 'rad_sw_out_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => rad_sw_out
+                ELSE
+                   to_be_resorted => rad_sw_out_av
+                ENDIF
 
              CASE ( 'rho_xy', 'rho_xz', 'rho_yz' )
Index: palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_3d.f90	(revision 1584)
+++ palm/trunk/SOURCE/data_output_3d.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Added support for RRTMG
 ! 
 ! Former revisions:
@@ -128,4 +128,9 @@
         
     USE pegrid
+
+    USE radiation_model_mod,                                                   &
+        ONLY:  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
+
 
     IMPLICIT NONE
@@ -467,4 +472,32 @@
              ELSE
                 to_be_resorted => qv_av
+             ENDIF
+
+          CASE ( 'rad_sw_in' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_sw_in
+             ELSE
+                to_be_resorted => rad_sw_in_av
+             ENDIF
+
+          CASE ( 'rad_sw_out' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_sw_out
+             ELSE
+                to_be_resorted => rad_sw_out_av
+             ENDIF
+
+          CASE ( 'rad_lw_in' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_lw_in
+             ELSE
+                to_be_resorted => rad_lw_in_av
+             ENDIF
+
+          CASE ( 'rad_lw_out' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_lw_out
+             ELSE
+                to_be_resorted => rad_lw_out_av
              ENDIF
 
Index: palm/trunk/SOURCE/data_output_mask.f90
===================================================================
--- palm/trunk/SOURCE/data_output_mask.f90	(revision 1584)
+++ palm/trunk/SOURCE/data_output_mask.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added support for RRTMG
 ! 
 ! Former revisions:
@@ -108,4 +108,8 @@
     
     USE pegrid
+
+    USE radiation_model_mod,                                                   &
+        ONLY:  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
 
     IMPLICIT NONE
@@ -393,4 +397,32 @@
              ENDIF
 
+          CASE ( 'rad_sw_in' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_sw_in
+             ELSE
+                to_be_resorted => rad_sw_in_av
+             ENDIF
+
+          CASE ( 'rad_sw_out' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_sw_out
+             ELSE
+                to_be_resorted => rad_sw_out_av
+             ENDIF
+
+          CASE ( 'rad_lw_in' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_lw_in
+             ELSE
+                to_be_resorted => rad_lw_in_av
+             ENDIF
+
+          CASE ( 'rad_lw_out' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_lw_out
+             ELSE
+                to_be_resorted => rad_lw_out_av
+             ENDIF
+
           CASE ( 'rho' )
              IF ( av == 0 )  THEN
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 1584)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 1585)
@@ -21,5 +21,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added output of timeseries and profiles for RRTMG
 ! 
 ! Former revisions:
@@ -178,6 +178,12 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY :  dots_rad, rad_net, rad_sw_in, radiation
-    
+        ONLY:  dots_rad, radiation, radiation_scheme, rad_net,                 &
+               rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out
+
+#if defined ( __rrtmg )
+    USE radiation_model_mod,                                                   &
+        ONLY:  rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
+#endif
+  
     USE statistics
 
@@ -741,5 +747,18 @@
              IF ( radiation )  THEN
                 sums_l(nzb,101,tn)  = sums_l(nzb,101,tn)  + rad_net(j,i)
-                sums_l(nzb,102,tn) = sums_l(nzb,102,tn) + rad_sw_in(j,i)
+                sums_l(nzb,102,tn)  = sums_l(nzb,102,tn)  + rad_lw_in(nzb,j,i)
+                sums_l(nzb,103,tn)  = sums_l(nzb,103,tn)  + rad_lw_out(nzb,j,i)
+                sums_l(nzb,104,tn)  = sums_l(nzb,104,tn)  + rad_sw_in(nzb,j,i)
+                sums_l(nzb,105,tn)  = sums_l(nzb,105,tn)  + rad_sw_out(nzb,j,i)
+
+#if defined ( __rrtmg )
+                IF ( radiation_scheme == 'rrtmg' )  THEN
+                   sums_l(nzb,106,tn)  = sums_l(nzb,106,tn)  + rrtm_aldif(0,j,i)
+                   sums_l(nzb,107,tn)  = sums_l(nzb,107,tn)  + rrtm_aldir(0,j,i)
+                   sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  + rrtm_asdif(0,j,i)
+                   sums_l(nzb,109,tn)  = sums_l(nzb,109,tn)  + rrtm_asdir(0,j,i)
+                ENDIF
+#endif
+
              ENDIF
 
@@ -1126,5 +1145,17 @@
        ENDIF
        
-
+       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
+          !$OMP DO
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                DO  k = nzb_s_inner(j,i)+1, nzt+1
+                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i) * rmask(j,i,sr)
+                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i) * rmask(j,i,sr)
+                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i) * rmask(j,i,sr)
+                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i) * rmask(j,i,sr)
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
 !
 !--    Calculate the user-defined profiles
@@ -1186,6 +1217,6 @@
           sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
           sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
-          sums(k,89:100)           = sums(k,89:100)     / ngp_2dh(sr)
-          sums(k,101:pr_palm-2)    = sums(k,101:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
+          sums(k,89:105)           = sums(k,89:105)     / ngp_2dh(sr)
+          sums(k,106:pr_palm-2)    = sums(k,106:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
        ENDDO
 
@@ -1319,6 +1350,19 @@
 
        IF ( radiation )  THEN
-          hom(:,1,101 ,sr) = sums(:,101)            ! rad_net
-          hom(:,1,102,sr)  = sums(:,102)            ! rad_sw_in
+          hom(:,1,101,sr) = sums(:,101)            ! rad_net
+          hom(:,1,102,sr) = sums(:,102)            ! rad_lw_in
+          hom(:,1,103,sr) = sums(:,103)            ! rad_lw_out
+          hom(:,1,104,sr) = sums(:,104)            ! rad_sw_in
+          hom(:,1,105,sr) = sums(:,105)            ! rad_sw_out
+
+#if defined ( __rrtmg )
+          IF ( radiation_scheme == 'rrtmg' )  THEN
+             hom(:,1,106,sr) = sums(:,106)            ! rrtm_aldif
+             hom(:,1,107,sr) = sums(:,107)            ! rrtm_aldir
+             hom(:,1,108,sr) = sums(:,108)            ! rrtm_asdif
+             hom(:,1,109,sr) = sums(:,109)            ! rrtm_asdir
+          ENDIF
+#endif
+
        ENDIF
 
@@ -1483,6 +1527,19 @@
 !--    Collect radiation model timeseries
        IF ( radiation )  THEN
-          ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)           ! rad_net
-          ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_sw_in
+          ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)          ! rad_net
+          ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_lw_in
+          ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr)          ! rad_lw_out
+          ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr)          ! rad_lw_in
+          ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_lw_out
+
+#if defined ( __rrtmg )
+          IF ( radiation_scheme == 'rrtmg' )  THEN
+             ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr)          ! rrtm_aldif
+             ts_value(dots_rad+6,sr) = hom(nzb,1,107,sr)          ! rrtm_aldir
+             ts_value(dots_rad+7,sr) = hom(nzb,1,108,sr)          ! rrtm_asdif
+             ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr)          ! rrtm_asdir
+          ENDIF
+#endif
+
        ENDIF
 
@@ -2951,4 +3008,18 @@
 
 
+       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
+          !$OMP DO
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                DO  k = nzb_s_inner(j,i)+1, nzt+1
+                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i) * rmask(j,i,sr)
+                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i) * rmask(j,i,sr)
+                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i) * rmask(j,i,sr)
+                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i) * rmask(j,i,sr)
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
+
 !
 !--    Calculate the user-defined profiles
@@ -3012,6 +3083,6 @@
           sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
           sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
-          sums(k,89:100)           = sums(k,89:100)     / ngp_2dh(sr)
-          sums(k,101:pr_palm-2)    = sums(k,101:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
+          sums(k,89:105)           = sums(k,89:105)     / ngp_2dh(sr)
+          sums(k,106:pr_palm-2)    = sums(k,106:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
        ENDDO
 
@@ -3288,6 +3359,19 @@
 !--    Collect radiation model timeseries
        IF ( radiation )  THEN
-          ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)           ! rad_net
-          ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_sw_in
+          ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)          ! rad_net
+          ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_lw_in
+          ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr)          ! rad_lw_out
+          ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr)          ! rad_lw_in
+          ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_lw_out
+
+#if defined ( __rrtmg )
+          IF ( radiation_scheme == 'rrtmg' )  THEN
+             ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr)          ! rrtm_aldif
+             ts_value(dots_rad+6,sr) = hom(nzb,1,107,sr)          ! rrtm_aldir
+             ts_value(dots_rad+7,sr) = hom(nzb,1,108,sr)          ! rrtm_asdif
+             ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr)          ! rrtm_asdir
+          ENDIF
+#endif
+
        ENDIF
 
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 1584)
+++ palm/trunk/SOURCE/header.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Further output for radiation model(s).
 !
 ! Former revisions:
@@ -245,6 +245,7 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  albedo, day_init, dt_radiation, lambda, net_radiation,          &
-               radiation, radiation_scheme, time_utc_init
+        ONLY:  albedo, albedo_type, albedo_type_name, constant_albedo,         &
+               day_init, dt_radiation, lambda, lw_radiation, net_radiation,    &
+               radiation, radiation_scheme, sw_radiation, time_utc_init
     
     USE spectrum,                                                              &
@@ -927,5 +928,5 @@
     IF ( radiation )  THEN
 !
-!--    Write land surface model header
+!--    Write radiation model header
        WRITE( io, 444 )
 
@@ -934,11 +935,19 @@
        ELSEIF ( radiation_scheme == "clear-sky" )  THEN
           WRITE( io, 446 )
-       ELSE
-          WRITE( io, 447 ) radiation_scheme
-       ENDIF
-
-       WRITE( io, 448 ) albedo
+       ELSEIF ( radiation_scheme == "rrtmg" )  THEN
+          WRITE( io, 447 )
+          IF ( .NOT. lw_radiation )  WRITE( io, 458 )
+          IF ( .NOT. sw_radiation )  WRITE( io, 459 )
+       ENDIF
+
+       IF ( albedo_type == 0 )  THEN
+          WRITE( io, 448 ) albedo
+       ELSE
+          WRITE( io, 456 ) albedo_type_name(albedo_type)
+       ENDIF
+       IF ( constant_albedo )  THEN
+          WRITE( io, 457 )
+       ENDIF
        WRITE( io, 449 ) dt_radiation
-
     ENDIF
 
@@ -2262,7 +2271,7 @@
 446 FORMAT (' --> Simple radiation scheme for clear sky is used (no clouds,',  &
                    ' default)')
-447 FORMAT (' --> Radiation scheme:', A)
-448 FORMAT (/'    Surface albedo: albedo = ', F5.3)
-449 FORMAT  ('    Timestep: dt_radiation = ', F5.2, '  s')
+447 FORMAT (' --> RRTMG scheme is used')
+448 FORMAT (/'     User-specific surface albedo: albedo = ', F5.3)
+449 FORMAT  ('     Timestep: dt_radiation = ', F5.2, '  s')
 
 450 FORMAT (//' LES / Turbulence quantities:'/ &
@@ -2273,4 +2282,8 @@
 454 FORMAT ('    TKE is not allowed to fall below ',E9.2,' (m/s)**2')
 455 FORMAT ('    initial TKE is prescribed as ',E9.2,' (m/s)**2')
+456 FORMAT (/'    Albedo is set for land surface type: ', A)
+457 FORMAT (/'    --> Albedo is fixed during the run')
+458 FORMAT (/'    --> Longwave radiation is disabled')
+459 FORMAT (/'    --> Shortwave radiation is disabled.')
 470 FORMAT (//' Actions during the simulation:'/ &
               ' -----------------------------'/)
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 1584)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Initialization of radiation code is now done after LSM initializtion
 !
 ! Former revisions:
@@ -1606,15 +1606,18 @@
     IF ( particle_advection )  CALL lpm_init
 
+!
+!-- If required, initialize quantities needed for the LSM
+    IF ( land_surface )  THEN
+       CALL location_message( 'initializing land surface model', .FALSE. )
+       CALL init_lsm
+       CALL location_message( 'finished', .TRUE. )
+    ENDIF
 
 !
 !-- If required, initialize radiation model
     IF ( radiation )  THEN
+       CALL location_message( 'initializing radiation model', .FALSE. )
        CALL init_radiation
-    ENDIF
-
-!
-!-- If required, initialize quantities needed for the LSM
-    IF ( land_surface )  THEN
-       CALL init_lsm
+       CALL location_message( 'finished', .TRUE. )
     ENDIF
 
Index: palm/trunk/SOURCE/land_surface_model.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model.f90	(revision 1584)
+++ palm/trunk/SOURCE/land_surface_model.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Modifications for RRTMG. Changed tables to PARAMETER type. 
 ! 
 ! Former revisions:
@@ -80,5 +80,5 @@
 !------------------------------------------------------------------------------!
      USE arrays_3d,                                                            &
-         ONLY:  pt, pt_p, q, q_p, qsws, rif, shf, ts, us, z0, z0h
+         ONLY:  pt, pt_p, q_p, qsws, rif, shf, ts, us, z0, z0h
 
      USE cloud_parameters,                                                     &
@@ -90,8 +90,8 @@
                 max_masks, precipitation, pt_surface, rho_surface,             &
                 roughness_length, surface_pressure, timestep_scheme, tsc,      &
-                z0h_factor
+                z0h_factor, time_since_reference_point
 
      USE indices,                                                              &
-         ONLY:  nxlg, nxrg, nyng, nysg, nzb_s_inner 
+         ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner 
 
      USE kinds
@@ -101,5 +101,5 @@
 
      USE radiation_model_mod,                                                  &
-         ONLY:  irad_scheme, rad_net, rad_sw_in, sigma_sb
+         ONLY:  radiation_scheme, rad_net, rad_sw_in, sigma_sb
 
      USE statistics,                                                           &
@@ -243,4 +243,5 @@
               qsws_veg_eb,      & !: surface flux of latent heat (vegetation portion)
               qsws_veg_eb_av,   & !: average of qsws_veg_eb
+              rad_net_l,        & !: local copy of rad_net (net radiation at surface)
               r_a,              & !: aerodynamic resistance 
               r_a_av,           & !: avergae of r_a
@@ -276,6 +277,4 @@
               t_soil_av, t_soil_1, t_soil_2,                                   &
               m_soil_av, m_soil_1, m_soil_2
-
-
 #endif
 
@@ -289,39 +288,39 @@
 !
 !-- Predefined Land surface classes (veg_type)
-    CHARACTER(26), DIMENSION(0:19) :: veg_type_name = (/          &
-                                   'user defined',                & ! 0 
-                                   'crops, mixed farming',        & !  1
-                                   'short grass',                 & !  2
-                                   'evergreen needleleaf trees',  & !  3
-                                   'deciduous needleleaf trees',  & !  4
-                                   'evergreen broadleaf trees' ,  & !  5
-                                   'deciduous broadleaf trees',   & !  6
-                                   'tall grass',                  & !  7
-                                   'desert',                      & !  8
-                                   'tundra',                      & !  9
-                                   'irrigated crops',             & ! 10
-                                   'semidesert',                  & ! 11
-                                   'ice caps and glaciers' ,      & ! 12
-                                   'bogs and marshes',            & ! 13
-                                   'inland water',                & ! 14
-                                   'ocean',                       & ! 15
-                                   'evergreen shrubs',            & ! 16
-                                   'deciduous shrubs',            & ! 17
-                                   'mixed forest/woodland',       & ! 18
-                                   'interrupted forest'           & ! 19
-                                                       /)
+    CHARACTER(26), DIMENSION(0:19), PARAMETER :: veg_type_name = (/ &
+                                   'user defined',                  & ! 0 
+                                   'crops, mixed farming',          & !  1
+                                   'short grass',                   & !  2
+                                   'evergreen needleleaf trees',    & !  3
+                                   'deciduous needleleaf trees',    & !  4
+                                   'evergreen broadleaf trees' ,    & !  5
+                                   'deciduous broadleaf trees',     & !  6
+                                   'tall grass',                    & !  7
+                                   'desert',                        & !  8
+                                   'tundra',                        & !  9
+                                   'irrigated crops',               & ! 10
+                                   'semidesert',                    & ! 11
+                                   'ice caps and glaciers' ,        & ! 12
+                                   'bogs and marshes',              & ! 13
+                                   'inland water',                  & ! 14
+                                   'ocean',                         & ! 15
+                                   'evergreen shrubs',              & ! 16
+                                   'deciduous shrubs',              & ! 17
+                                   'mixed forest/woodland',         & ! 18
+                                   'interrupted forest'             & ! 19
+                                                                 /)
 
 !
 !-- Soil model classes (soil_type)
-    CHARACTER(12), DIMENSION(0:7) :: soil_type_name = (/ &
-                                   'user defined',        & ! 0 
-                                   'coarse',              & ! 1
-                                   'medium',              & ! 2
-                                   'medium-fine',         & ! 3
-                                   'fine',                & ! 4
-                                   'very fine' ,          & ! 5
-                                   'organic',             & ! 6
-                                   'loamy (CH)'           & ! 7
-                                                        /)
+    CHARACTER(12), DIMENSION(0:7), PARAMETER :: soil_type_name = (/ &
+                                   'user defined',                  & ! 0 
+                                   'coarse',                        & ! 1
+                                   'medium',                        & ! 2
+                                   'medium-fine',                   & ! 3
+                                   'fine',                          & ! 4
+                                   'very fine' ,                    & ! 5
+                                   'organic',                       & ! 6
+                                   'loamy (CH)'                     & ! 7
+                                                                 /)
 !
 !-- Land surface parameters according to the respective classes (veg_type)
@@ -330,96 +329,96 @@
 !-- Land surface parameters I
 !--                          r_canopy_min,     lai,   c_veg,     g_d
-    REAL(wp), DIMENSION(0:3,1:19) :: veg_pars = RESHAPE( (/           &
-                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & !  1
-                                 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp, & !  2
-                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & !  3
-                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & !  4
-                                 175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & !  5
-                                 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp, & !  6
-                                 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp, & !  7
-                                 250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp, & !  8
-                                  80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp, & !  9
-                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 10
-                                 150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp, & ! 11
-                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12
-                                 240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp, & ! 13
-                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14
-                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15
-                                 225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp, & ! 16
-                                 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp, & ! 17
-                                 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 18
-                                 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp  & ! 19
+    REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: veg_pars = RESHAPE( (/ &
+                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & !  1
+                                 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,  & !  2
+                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  3
+                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  4
+                                 175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  5
+                                 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,  & !  6
+                                 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,  & !  7
+                                 250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  & !  8
+                                  80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  & !  9
+                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & ! 10
+                                 150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,  & ! 11
+                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 12
+                                 240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,  & ! 13
+                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 14
+                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 15
+                                 225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,  & ! 16
+                                 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,  & ! 17
+                                 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & ! 18
+                                 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp   & ! 19
                                  /), (/ 4, 19 /) )
 
 !
 !-- Land surface parameters II          z0,         z0h
-    REAL(wp), DIMENSION(0:1,1:19) :: roughness_par = RESHAPE( (/ & 
-                                   0.25_wp,  0.25E-2_wp,         & !  1
-                                   0.20_wp,  0.20E-2_wp,         & !  2
-                                   2.00_wp,     2.00_wp,         & !  3
-                                   2.00_wp,     2.00_wp,         & !  4
-                                   2.00_wp,     2.00_wp,         & !  5
-                                   2.00_wp,     2.00_wp,         & !  6
-                                   0.47_wp,  0.47E-2_wp,         & !  7
-                                  0.013_wp, 0.013E-2_wp,         & !  8
-                                  0.034_wp, 0.034E-2_wp,         & !  9
-                                    0.5_wp,  0.50E-2_wp,         & ! 10
-                                   0.17_wp,  0.17E-2_wp,         & ! 11
-                                 1.3E-3_wp,   1.3E-4_wp,         & ! 12
-                                   0.83_wp,  0.83E-2_wp,         & ! 13
-                                   0.00_wp,  0.00E-2_wp,         & ! 14
-                                   0.00_wp,  0.00E-2_wp,         & ! 15
-                                   0.10_wp,  0.10E-2_wp,         & ! 16
-                                   0.25_wp,  0.25E-2_wp,         & ! 17
-                                   2.00_wp,  2.00E-2_wp,         & ! 18
-                                   1.10_wp,  1.10E-2_wp          & ! 19
+    REAL(wp), DIMENSION(0:1,1:19), PARAMETER :: roughness_par = RESHAPE( (/ & 
+                                   0.25_wp,  0.25E-2_wp,                    & !  1
+                                   0.20_wp,  0.20E-2_wp,                    & !  2
+                                   2.00_wp,     2.00_wp,                    & !  3
+                                   2.00_wp,     2.00_wp,                    & !  4
+                                   2.00_wp,     2.00_wp,                    & !  5
+                                   2.00_wp,     2.00_wp,                    & !  6
+                                   0.47_wp,  0.47E-2_wp,                    & !  7
+                                  0.013_wp, 0.013E-2_wp,                    & !  8
+                                  0.034_wp, 0.034E-2_wp,                    & !  9
+                                    0.5_wp,  0.50E-2_wp,                    & ! 10
+                                   0.17_wp,  0.17E-2_wp,                    & ! 11
+                                 1.3E-3_wp,   1.3E-4_wp,                    & ! 12
+                                   0.83_wp,  0.83E-2_wp,                    & ! 13
+                                   0.00_wp,  0.00E-2_wp,                    & ! 14
+                                   0.00_wp,  0.00E-2_wp,                    & ! 15
+                                   0.10_wp,  0.10E-2_wp,                    & ! 16
+                                   0.25_wp,  0.25E-2_wp,                    & ! 17
+                                   2.00_wp,  2.00E-2_wp,                    & ! 18
+                                   1.10_wp,  1.10E-2_wp                     & ! 19
                                  /), (/ 2, 19 /) )
 
 !
 !-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
-    REAL(wp), DIMENSION(0:2,1:19) :: surface_pars = RESHAPE( (/           &
-                                      10.0_wp,       10.0_wp, 0.05_wp, & !  1
-                                      10.0_wp,       10.0_wp, 0.05_wp, & !  2
-                                      20.0_wp,       15.0_wp, 0.03_wp, & !  3
-                                      20.0_wp,       15.0_wp, 0.03_wp, & !  4
-                                      20.0_wp,       15.0_wp, 0.03_wp, & !  5
-                                      20.0_wp,       15.0_wp, 0.03_wp, & !  6
-                                      10.0_wp,       10.0_wp, 0.05_wp, & !  7
-                                      15.0_wp,       15.0_wp, 0.00_wp, & !  8
-                                      10.0_wp,       10.0_wp, 0.05_wp, & !  9
-                                      10.0_wp,       10.0_wp, 0.05_wp, & ! 10 
-                                      10.0_wp,       10.0_wp, 0.05_wp, & ! 11
-                                      58.0_wp,       58.0_wp, 0.00_wp, & ! 12
-                                      10.0_wp,       10.0_wp, 0.05_wp, & ! 13 
-                                    1.0E20_wp,     1.0E20_wp, 0.00_wp, & ! 14
-                                    1.0E20_wp,     1.0E20_wp, 0.00_wp, & ! 15
-                                      10.0_wp,       10.0_wp, 0.05_wp, & ! 16
-                                      10.0_wp,       10.0_wp, 0.05_wp, & ! 17
-                                      20.0_wp,       15.0_wp, 0.03_wp, & ! 18
-                                      20.0_wp,       15.0_wp, 0.03_wp  & ! 19
+    REAL(wp), DIMENSION(0:2,1:19), PARAMETER :: surface_pars = RESHAPE( (/ &
+                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  1
+                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  2
+                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  3
+                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  4
+                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  5
+                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  6
+                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  7
+                                      15.0_wp,       15.0_wp, 0.00_wp,     & !  8
+                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  9
+                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 10 
+                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 11
+                                      58.0_wp,       58.0_wp, 0.00_wp,     & ! 12
+                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 13 
+                                    1.0E20_wp,     1.0E20_wp, 0.00_wp,     & ! 14
+                                    1.0E20_wp,     1.0E20_wp, 0.00_wp,     & ! 15
+                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 16
+                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 17
+                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 18
+                                      20.0_wp,       15.0_wp, 0.03_wp      & ! 19
                                       /), (/ 3, 19 /) )
 
 !
 !-- Root distribution (sum = 1)  level 1, level 2, level 3, level 4,
-    REAL(wp), DIMENSION(0:3,1:19) :: root_distribution = RESHAPE( (/ &
-                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & !  1
-                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & !  2
-                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp, & !  3
-                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp, & !  4
-                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp, & !  5
-                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp, & !  6
-                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp, & !  7
-                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & !  8
-                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp, & !  9
-                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 10
-                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp, & ! 11
-                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12
-                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp, & ! 13
-                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14
-                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15
-                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 16 
-                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 17 
-                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 18
-                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp  & ! 19
+    REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: root_distribution = RESHAPE( (/ &
+                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  1
+                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  2
+                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp,            & !  3
+                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp,            & !  4
+                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp,            & !  5
+                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp,            & !  6
+                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp,            & !  7
+                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  8
+                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp,            & !  9
+                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & ! 10
+                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp,            & ! 11
+                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 12
+                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp,            & ! 13
+                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 14
+                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 15
+                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16 
+                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 17 
+                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 18
+                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp             & ! 19
                                  /), (/ 4, 19 /) )
 
@@ -429,5 +428,5 @@
 !
 !-- Soil parameters I           alpha_vg,      l_vg,    n_vg, gamma_w_sat
-    REAL(wp), DIMENSION(0:3,1:7) :: soil_pars = RESHAPE( (/                &
+    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/     &
                                  3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, & ! 1
                                  3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, & ! 2
@@ -441,5 +440,5 @@
 !
 !-- Soil parameters II              m_sat,     m_fc,   m_wilt,    m_res  
-    REAL(wp), DIMENSION(0:3,1:7) :: m_soil_pars = RESHAPE( (/            &
+    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
                                  0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
                                  0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
@@ -490,5 +489,4 @@
 !-- Public prognostic variables
     PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
-
 
     INTERFACE init_lsm
@@ -573,4 +571,5 @@
        ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) )
        ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( rad_net_l(nysg:nyng,nxlg:nxrg) )
        ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
        ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
@@ -583,5 +582,5 @@
 #if ! defined( __nopointer )
 !
-!-- Initial assignment of the pointers
+!--    Initial assignment of the pointers
        t_soil    => t_soil_1;    t_soil_p    => t_soil_2
        t_surface => t_surface_1; t_surface_p => t_surface_2
@@ -890,4 +889,5 @@
 !
 !--    Add timeseries for land surface model
+
        dots_label(dots_num+1) = "ghf_eb"
        dots_label(dots_num+2) = "shf_eb"
@@ -904,5 +904,4 @@
        dots_soil = dots_num + 1
        dots_num  = dots_num + 8
-
 
        RETURN
@@ -957,4 +956,5 @@
        exn = ( surface_pressure / 1000.0_wp )**0.286_wp
 
+       
 
        DO i = nxlg, nxrg
@@ -984,7 +984,8 @@
 !--          f1: correction for incoming shortwave radiation (stomata close at 
 !--          night)
-             IF ( irad_scheme /= 0 )  THEN
-                f1 = MIN(1.0_wp, ( 0.004_wp * rad_sw_in(j,i) + 0.05_wp ) /     &
-                              (0.81_wp * (0.004_wp * rad_sw_in(j,i) + 1.0_wp) ))
+             IF ( radiation_scheme /= 'constant' )  THEN
+                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /&
+                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)        &
+                               + 1.0_wp)) )
              ELSE
                 f1 = 1.0_wp
@@ -1112,5 +1113,5 @@
 !
 !--          Add LW up so that it can be removed in prognostic equation
-             rad_net(j,i) = rad_net(j,i) + sigma_sb * t_surface(j,i) ** 4
+             rad_net_l(j,i) = rad_net(j,i) + sigma_sb * t_surface(j,i) ** 4
 
              IF ( humidity )  THEN
@@ -1118,5 +1119,5 @@
 !
 !--             Numerator of the prognostic equation
-                coef_1 = rad_net(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
+                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
                          + f_shf / exn * T_1 + f_qsws * ( q_p(k+1,j,i) - q_s   &
                          + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
@@ -1132,5 +1133,5 @@
 !
 !--             Numerator of the prognostic equation
-                coef_1 = rad_net(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
+                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
                          + f_shf / exn * T_1 + lambda_surface                  &
                          * t_soil(nzb_soil,j,i)
@@ -1175,5 +1176,5 @@
 !
 !--          Calculate fluxes
-             rad_net(j,i)   = rad_net(j,i) + 3.0_wp * sigma_sb                 &
+             rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb                 &
                               * t_surface(j,i)**4 - 4.0_wp * sigma_sb          &
                               * t_surface(j,i)**3 * t_surface_p(j,i)
@@ -1621,14 +1622,4 @@
           CASE ( 0 )
 
-             t_surface = t_surface_p
-             t_soil    = t_soil_p
-             IF ( humidity )  THEN
-                m_soil    = m_soil_p
-                m_liq_eb  = m_liq_eb_p
-             ENDIF
-
-
-          CASE ( 1 )
-
              t_surface  => t_surface_1; t_surface_p  => t_surface_2
              t_soil     => t_soil_1;    t_soil_p     => t_soil_2
@@ -1638,4 +1629,14 @@
              ENDIF
 
+
+          CASE ( 1 )
+
+             t_surface  => t_surface_2; t_surface_p  => t_surface_1
+             t_soil     => t_soil_2;    t_soil_p     => t_soil_1
+             IF ( humidity )  THEN
+                m_soil    => m_soil_2;   m_soil_p    => m_soil_1
+                m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
+             ENDIF
+
        END SELECT
 #endif
Index: palm/trunk/SOURCE/netcdf.f90
===================================================================
--- palm/trunk/SOURCE/netcdf.f90	(revision 1584)
+++ palm/trunk/SOURCE/netcdf.f90	(revision 1585)
@@ -148,4 +148,7 @@
     USE profil_parameter,                                                       &
         ONLY:  crmax, cross_profiles, dopr_index,profile_columns, profile_rows
+
+    USE radiation_model_mod,                                                    &
+        ONLY:  rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out
 
     USE spectrum,                                                               &
@@ -523,5 +526,6 @@
                 CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',&
                        'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',    &
-                       'rho', 's', 'sa', 'vpt' )
+                       'rho', 's', 'sa', 'vpt', 'rad_lw_in', 'rad_lw_out', &
+                       'rad_sw_in', 'rad_sw_out' )
 
                    grid_x = 'x'
@@ -1103,5 +1107,6 @@
                 CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',   &
                        'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho',&
-                       's', 'sa', 'vpt' )
+                       's', 'sa', 'vpt' , 'rad_lw_in', 'rad_lw_out',          &
+                       'rad_sw_in', 'rad_sw_out' )
 
                    grid_x = 'x'
@@ -1756,5 +1761,6 @@
                              'prr_xy', 'pt_xy', 'q_xy', 'qc_xy', 'ql_xy',        &
                              'ql_c_xy', 'ql_v_xy', 'ql_vp_xy', 'qr_xy', 'qv_xy', &
-                             'rho_xy', 's_xy', 'sa_xy', 'vpt_xy' )
+                             'rho_xy', 's_xy', 'sa_xy', 'vpt_xy', 'rad_lw_in_xy',&
+                             'rad_lw_out_xy', 'rad_sw_in_xy', 'rad_sw_out_xy' )
 
                          grid_x = 'x'
@@ -2439,5 +2445,6 @@
                           'prr_xz', 'pt_xz', 'q_xz', 'qc_xz', 'ql_xz',         &
                           'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz', 'qv_xz',  &
-                          'rho_xz', 's_xz', 'sa_xz', 'vpt_xz' )
+                          'rho_xz', 's_xz', 'sa_xz', 'vpt_xz' , 'rad_lw_in_xz',&
+                             'rad_lw_out_xz', 'rad_sw_in_xz', 'rad_sw_out_xz' )
 
                       grid_x = 'x'
@@ -3115,5 +3122,6 @@
                           'prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz',        &
                           'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz', &
-                          'rho_yz', 's_yz', 'sa_yz', 'vpt_yz' )
+                          'rho_yz', 's_yz', 'sa_yz', 'vpt_yz', 'rad_lw_in_yz',&
+                             'rad_lw_out_yz', 'rad_sw_in_yz', 'rad_sw_out_yz' )
 
                       grid_x = 'x'
Index: palm/trunk/SOURCE/package_parin.f90
===================================================================
--- palm/trunk/SOURCE/package_parin.f90	(revision 1584)
+++ palm/trunk/SOURCE/package_parin.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added several radiation_par parameters
 !
 ! Former revisions:
@@ -144,6 +144,8 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY: albedo, day_init, dt_radiation, lambda, net_radiation, radiation,&
-              radiation_scheme, time_radiation, time_utc_init 
+        ONLY: albedo, albedo_type, albedo_lw_dif, albedo_lw_dir, albedo_sw_dif,&
+              albedo_sw_dir, day_init, constant_albedo, dt_radiation, lambda,  &
+              lw_radiation, net_radiation, radiation, radiation_scheme,        &
+              sw_radiation, time_utc_init
                 
 
@@ -221,6 +223,9 @@
                                   write_particle_statistics
 
-    NAMELIST /radiation_par/      albedo, day_init, dt_radiation, lambda,      &
-                                  net_radiation, radiation_scheme, time_utc_init 
+    NAMELIST /radiation_par/      albedo, albedo_type, albedo_lw_dir,          &
+                                  albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
+                                  constant_albedo, day_init, dt_radiation,     &
+                                  lambda, lw_radiation, net_radiation,         &
+                                  radiation_scheme, sw_radiation, time_utc_init 
 
     NAMELIST /spectra_par/        averaging_interval_sp, comp_spectra_level,   &
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 1584)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! Added call for temperature tendency calculation due to radiative flux divergence
 ! 
 ! Former revisions:
@@ -277,4 +277,7 @@
         ONLY:  production_e, production_e_acc
 
+    USE radiation_model_mod,                                                   &
+        ONLY:  radiation, radiation_scheme, radiation_tendency
+
     USE statistics,                                                            &
         ONLY:  hom
@@ -593,4 +596,11 @@
                 CALL subsidence( i, j, tend, pt, pt_init, 2 )
              ENDIF
+
+!
+!--          If required, add tendency due to radiative heating/cooling
+             IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
+                CALL radiation_tendency ( i, j, tend )
+             ENDIF
+
 
              CALL user_actions( i, j, 'pt-tendency' )
@@ -1246,4 +1256,10 @@
             .NOT. use_subsidence_tendencies )  THEN
           CALL subsidence( tend, pt, pt_init, 2 )
+       ENDIF
+
+!
+!--    If required, add tendency due to radiative heating/cooling
+       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
+          CALL radiation_tendency ( tend )
        ENDIF
 
Index: palm/trunk/SOURCE/radiation_model.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model.f90	(revision 1584)
+++ 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
Index: palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/read_3d_binary.f90	(revision 1584)
+++ palm/trunk/SOURCE/read_3d_binary.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adapted for RRTMG
 ! 
 ! Former revisions:
@@ -111,5 +111,6 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY: rad_net_av, rad_sw_in_av
+        ONLY: rad_net, rad_net_av, 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
 
     USE random_function_mod,                                                   &
@@ -333,5 +334,5 @@
 !--    First compare the version numbers
        READ ( 13 )  version_on_file
-       binary_version = '4.0'
+       binary_version = '4.1'
        IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
           WRITE( message_string, * ) 'version mismatch concerning data ',      &
@@ -812,4 +813,12 @@
                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
+                CASE ( 'rad_net' )
+                   IF ( .NOT. ALLOCATED( rad_net ) )  THEN
+                      ALLOCATE( rad_net(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
                 CASE ( 'rad_net_av' )
                    IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
@@ -819,12 +828,67 @@
                    rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                CASE ( 'rad_lw_in' )
+                   IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
+                      ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_in_av' )
+                   IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
+                      ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_out' )
+                   IF ( .NOT. ALLOCATED( rad_lw_out ) )  THEN
+                      ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_out_av' )
+                   IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
+                      ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_in' )
+                   IF ( .NOT. ALLOCATED( rad_sw_in ) )  THEN
+                      ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'rad_sw_in_av' )
                    IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
-                      ALLOCATE( rad_sw_in_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF  
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   rad_sw_in_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
-                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                      ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_out' )
+                   IF ( .NOT. ALLOCATED( rad_sw_out ) )  THEN
+                      ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_out_av' )
+                   IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
+                      ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'random_iv' )  ! still unresolved issue
Index: palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- palm/trunk/SOURCE/read_var_list.f90	(revision 1584)
+++ palm/trunk/SOURCE/read_var_list.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Adapted for RRTMG
 ! 
 ! Former revisions:
@@ -162,4 +162,7 @@
 
     USE pegrid
+
+    USE radiation_model_mod,                                                   &
+        ONLY: time_radiation
 
     USE statistics,                                                            &
@@ -609,4 +612,6 @@
           CASE ( 'time_dvrp' )
              READ ( 13 )  time_dvrp
+          CASE ( 'time_radiation' )
+             READ ( 13 )  time_radiation
           CASE ( 'time_restart' )
              READ ( 13 )  time_restart
Index: palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 1584)
+++ palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adapted for RRTMG
 ! 
 ! Former revisions:
@@ -112,5 +112,7 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av
+        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av, rad_sw_out,       &
+               rad_sw_out_av, rad_lw_in, rad_lw_in_av, rad_lw_out,             &
+               rad_lw_out_av
 
     IMPLICIT NONE
@@ -326,9 +328,27 @@
                 rad_net_av = 0.0_wp
 
-             CASE ( 'rad_sw_in*' )
+             CASE ( 'rad_lw_in' )
+                IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
+                   ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_lw_in_av = 0.0_wp
+
+             CASE ( 'rad_lw_out' )
+                IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
+                   ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_lw_out_av = 0.0_wp
+
+             CASE ( 'rad_sw_in' )
                 IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
-                   ALLOCATE( rad_sw_in_av(nysg:nyng,nxlg:nxrg) )
+                   ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 ENDIF
                 rad_sw_in_av = 0.0_wp
+
+             CASE ( 'rad_sw_out' )
+                IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
+                   ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_sw_out_av = 0.0_wp
 
              CASE ( 'rho' )
@@ -732,8 +752,38 @@
              ENDDO
 
-          CASE ( 'rad_sw_in*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   rad_sw_in_av(j,i) = rad_sw_in_av(j,i) + rad_sw_in(j,i)
+          CASE ( 'rad_lw_in' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) + rad_lw_in(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_lw_out' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) + rad_lw_out(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+
+          CASE ( 'rad_sw_in' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) + rad_sw_in(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_out' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) + rad_sw_out(k,j,i)
+                   ENDDO
                 ENDDO
              ENDDO
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 1584)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Moved call of radiation scheme. Added support for RRTM
 ! 
 ! Former revisions:
@@ -226,6 +226,6 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY: dt_radiation, irad_scheme, radiation, radiation_clearsky,        &
-              radiation_constant, radiation_rrtm, time_radiation
+        ONLY: dt_radiation, radiation, radiation_clearsky,                     &
+              radiation_rrtmg, radiation_scheme, time_radiation
 
     USE statistics,                                                            &
@@ -362,4 +362,23 @@
           ENDIF
 
+          IF ( radiation .AND. intermediate_timestep_count                     &
+               == intermediate_timestep_count_max  )  THEN
+
+               time_radiation = time_radiation + dt_3d
+
+             IF ( time_radiation >= dt_radiation )  THEN
+
+                CALL cpu_log( log_point(50), 'radiation', 'start' )
+
+                time_radiation = time_radiation - dt_radiation
+                IF ( radiation_scheme == 'clear-sky' )  THEN
+                   CALL radiation_clearsky
+                ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
+                   CALL radiation_rrtmg
+                ENDIF
+
+                CALL cpu_log( log_point(50), 'radiation', 'stop' )
+             ENDIF
+          ENDIF
 !
 !--       Solve the prognostic equations. A fast cache optimized version with
@@ -597,22 +616,4 @@
 !--       2) run soil model 
           IF ( land_surface )  THEN
-
-             time_radiation = time_radiation + dt_3d
-
-             IF ( time_radiation >= dt_radiation )  THEN
-
-                CALL cpu_log( log_point(50), 'radiation', 'start' )
-
-                time_radiation = time_radiation - dt_radiation
-                IF ( irad_scheme == 0 )  THEN
-                   CALL radiation_constant              
-                ELSEIF ( irad_scheme == 1 )  THEN
-                   CALL radiation_clearsky
-                ELSEIF ( irad_scheme == 2 )  THEN
-                   CALL radiation_rrtm
-                ENDIF
-
-                CALL cpu_log( log_point(50), 'radiation', 'stop' )
-             ENDIF
 
              CALL cpu_log( log_point(54), 'land_surface', 'start' )
@@ -693,4 +694,5 @@
              !$acc update device( vpt ) 
           ENDIF 
+
 !
 !--       Compute the diffusion quantities
Index: palm/trunk/SOURCE/user_init_land_surface.f90
===================================================================
--- palm/trunk/SOURCE/user_init_land_surface.f90	(revision 1584)
+++ palm/trunk/SOURCE/user_init_land_surface.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Changed description text
 ! 
 ! Former revisions:
@@ -52,5 +52,5 @@
 
 !
-!-- Here the user-defined grid initializing actions follow:
+!-- Here the user-defined land surface initialization actions follow:
 
 
Index: palm/trunk/SOURCE/user_init_radiation.f90
===================================================================
--- palm/trunk/SOURCE/user_init_radiation.f90	(revision 1585)
+++ palm/trunk/SOURCE/user_init_radiation.f90	(revision 1585)
@@ -0,0 +1,56 @@
+ SUBROUTINE user_init_radiation
+
+!--------------------------------------------------------------------------------!
+! This file is part of PALM.
+!
+! 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 1997-2014 Leibniz Universitaet Hannover
+!--------------------------------------------------------------------------------!
+!
+! Current revisions:
+! -----------------
+! Initial revision
+! 
+! Former revisions:
+! -----------------
+! $Id$
+! 
+! Description:
+! ------------
+! Execution of user-defined actions to initiate the radiation model
+!------------------------------------------------------------------------------!
+
+    USE arrays_3d
+    
+    USE control_parameters
+    
+    USE indices
+    
+    USE kinds
+
+    USE radiation_model_mod
+    
+    USE user
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) :: i   !: running index
+    INTEGER(iwp) :: j   !: running index
+
+!
+!-- Here the user-defined radiation initialization actions follow:
+
+
+
+ END SUBROUTINE user_init_radiation
+
Index: palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/write_3d_binary.f90	(revision 1584)
+++ palm/trunk/SOURCE/write_3d_binary.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adapted for RRTMG
 ! 
 ! Former revisions:
@@ -102,5 +102,7 @@
     
     USE radiation_model_mod,                                                   &
-        ONLY: radiation, rad_net_av, rad_sw_in_av
+        ONLY: radiation, rad_net, rad_net_av, 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
 
     USE random_function_mod,                                                   &
@@ -124,5 +126,5 @@
 !
 !-- Write arrays.
-    binary_version = '4.0'
+    binary_version = '4.1'
 
     WRITE ( 14 )  binary_version
@@ -262,10 +264,34 @@
        ENDIF 
     ENDIF
+    IF ( ALLOCATED( rad_net ) )  THEN
+       WRITE ( 14 )  'rad_net             ';  WRITE ( 14 )  rad_net  
+    ENDIF
     IF ( radiation )  THEN
        IF ( ALLOCATED( rad_net_av ) )  THEN
           WRITE ( 14 )  'rad_net_av          ';  WRITE ( 14 )  rad_net_av  
        ENDIF  
+       IF ( ALLOCATED( rad_lw_in ) )  THEN
+          WRITE ( 14 )  'rad_lw_in           ';  WRITE ( 14 )  rad_lw_in  
+       ENDIF
+       IF ( ALLOCATED( rad_lw_in_av ) )  THEN
+          WRITE ( 14 )  'rad_lw_in_av        ';  WRITE ( 14 )  rad_lw_in_av  
+       ENDIF 
+       IF ( ALLOCATED( rad_lw_out ) )  THEN
+          WRITE ( 14 )  'rad_lw_out          ';  WRITE ( 14 )  rad_lw_out 
+       ENDIF 
+       IF ( ALLOCATED( rad_lw_out_av ) )  THEN
+          WRITE ( 14 )  'rad_lw_out_av       ';  WRITE ( 14 )  rad_lw_out_av  
+       ENDIF 
+       IF ( ALLOCATED( rad_sw_in ) )  THEN
+          WRITE ( 14 )  'rad_sw_in           ';  WRITE ( 14 )  rad_sw_in  
+       ENDIF 
        IF ( ALLOCATED( rad_sw_in_av ) )  THEN
-          WRITE ( 14 )  'rad_sw_in_av          ';  WRITE ( 14 )  rad_sw_in_av  
+          WRITE ( 14 )  'rad_sw_in_av        ';  WRITE ( 14 )  rad_sw_in_av  
+       ENDIF 
+       IF ( ALLOCATED( rad_sw_out ) )  THEN
+          WRITE ( 14 )  'rad_sw_out          ';  WRITE ( 14 )  rad_sw_out  
+       ENDIF 
+       IF ( ALLOCATED( rad_sw_out_av ) )  THEN
+          WRITE ( 14 )  'rad_sw_out_av       ';  WRITE ( 14 )  rad_sw_out_av  
        ENDIF 
     ENDIF
Index: palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- palm/trunk/SOURCE/write_var_list.f90	(revision 1584)
+++ palm/trunk/SOURCE/write_var_list.f90	(revision 1585)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adapted for RRTMG
 ! 
 ! Former revisions:
@@ -147,4 +147,7 @@
     
     USE pegrid
+
+    USE radiation_model_mod,                                                   &
+        ONLY: time_radiation
 
     USE statistics,                                                            &
@@ -522,4 +525,6 @@
     WRITE ( 14 )  'time_dvrp                     '
     WRITE ( 14 )  time_dvrp
+    WRITE ( 14 )  'time_radiation                '
+    WRITE ( 14 )  time_radiation
     WRITE ( 14 )  'time_restart                  '
     WRITE ( 14 )  time_restart
