Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 1975)
+++ palm/trunk/SOURCE/Makefile	(revision 1976)
@@ -509,5 +509,5 @@
 mod_particle_attributes.o: mod_particle_attributes.f90 mod_kinds.o
 netcdf_interface_mod.o: netcdf_interface_mod.f90 modules.o mod_kinds.o \
-   land_surface_model_mod.o spectra_mod.o
+   land_surface_model_mod.o radiation_model_mod.o spectra_mod.o
 nudging_mod.o: modules.o cpulog_mod.o mod_kinds.o
 package_parin.o: modules.o mod_kinds.o mod_particle_attributes.o
Index: palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/average_3d_data.f90	(revision 1975)
+++ palm/trunk/SOURCE/average_3d_data.f90	(revision 1976)
@@ -101,9 +101,5 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  rad_net, rad_net_av, rad_lw_in, rad_lw_in_av, rad_lw_out,       &
-               rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr,        &
-               rad_lw_hr_av, rad_sw_in, rad_sw_in_av, rad_sw_out,              &
-               rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,        &
-               rad_sw_hr_av
+        ONLY:  radiation, radiation_3d_data_averaging
 
 
@@ -295,83 +291,4 @@
              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_lw_cs_hr' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
-                   ENDDO
-                ENDDO
-             ENDDO
-
-          CASE ( 'rad_lw_hr' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      rad_lw_hr_av(k,j,i) = rad_lw_hr_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
-
-          CASE ( 'rad_sw_cs_hr' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
-                   ENDDO
-                ENDDO
-             ENDDO
-
-          CASE ( 'rad_sw_hr' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
-                   ENDDO
-                ENDDO
-             ENDDO
-
           CASE ( 'rho' )
              DO  i = nxlg, nxrg
@@ -487,4 +404,10 @@
 
 !
+!--          Radiation quantity
+             IF ( radiation )  THEN
+                CALL radiation_3d_data_averaging( 'average', doav(ii) )
+             ENDIF
+
+!
 !--          User-defined quantity
              CALL user_3d_data_averaging( 'average', doav(ii) )
Index: palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 1975)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 1976)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Output of radiation quantities is now done directly in the respective module
 ! 
 ! Former revisions:
@@ -26,5 +26,6 @@
 !
 ! 1972 2016-07-26 07:52:02Z maronga
-! Output of land surface quantities is now done directly in the respective module
+! Output of land surface quantities is now done directly in the respective 
+! module
 ! 
 ! 1960 2016-07-12 16:34:24Z suehring
@@ -195,9 +196,5 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av, rad_sw_out,       &
-               rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,        &
-               rad_sw_hr_av, rad_lw_in, rad_lw_in_av, rad_lw_out,              &
-               rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr,        &
-               rad_lw_hr_av
+        ONLY:  radiation, radiation_data_output_2d
 
     IMPLICIT NONE
@@ -753,86 +750,5 @@
                 IF ( mode == 'xy' )  level_z = zu
 
-             CASE ( 'rad_net*_xy' )        ! 2d-array
-                IF ( av == 0 ) THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) =  rad_net(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng 
-                         local_pf(i,j,nzb+1) =  rad_net_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
-                IF ( mode == 'xy' )  level_z = zu
-
-             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
-                IF ( mode == 'xy' )  level_z = zu
-
-             CASE ( 'rad_lw_cs_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_cs_hr_yz' )
-                IF ( av == 0 )  THEN
-                   to_be_resorted => rad_lw_cs_hr
-                ELSE
-                   to_be_resorted => rad_lw_cs_hr_av
-                ENDIF
-                IF ( mode == 'xy' )  level_z = zw
-
-             CASE ( 'rad_lw_hr_xy', 'rad_lw_hr_xz', 'rad_lw_hr_yz' )
-                IF ( av == 0 )  THEN
-                   to_be_resorted => rad_lw_hr
-                ELSE
-                   to_be_resorted => rad_lw_hr_av
-                ENDIF
-                IF ( mode == 'xy' )  level_z = zw
-
-             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
-                IF ( mode == 'xy' )  level_z = zu
-
-             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
-                IF ( mode == 'xy' )  level_z = zu
-
-             CASE ( 'rad_sw_cs_hr_xy', 'rad_sw_cs_hr_xz', 'rad_sw_cs_hr_yz' )
-                IF ( av == 0 )  THEN
-                   to_be_resorted => rad_sw_cs_hr
-                ELSE
-                   to_be_resorted => rad_sw_cs_hr_av
-                ENDIF
-                IF ( mode == 'xy' )  level_z = zw
-
-             CASE ( 'rad_sw_hr_xy', 'rad_sw_hr_xz', 'rad_sw_hr_yz' )
-                IF ( av == 0 )  THEN
-                   to_be_resorted => rad_sw_hr
-                ELSE
-                   to_be_resorted => rad_sw_hr_av
-                ENDIF
-                IF ( mode == 'xy' )  level_z = zw
+
 
              CASE ( 'rho_xy', 'rho_xz', 'rho_yz' )
@@ -1034,4 +950,11 @@
                    CALL lsm_data_output_2d( av, do2d(av,if), found, grid, mode,&
                                             local_pf, two_d, nzb_do, nzt_do )
+                ENDIF
+
+!
+!--             Radiation quantity
+                IF ( .NOT. found  .AND.  radiation )  THEN
+                   CALL radiation_data_output_2d( av, do2d(av,if), found, grid,&
+                                                  mode, local_pf, two_d  )
                 ENDIF
 
Index: palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_3d.f90	(revision 1975)
+++ palm/trunk/SOURCE/data_output_3d.f90	(revision 1976)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Output of radiation quantities is now done directly in the respective module
 ! 
 ! Former revisions:
@@ -165,8 +165,5 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  rad_lw_in, rad_lw_in_av, rad_lw_out, rad_lw_out_av,             &
-               rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av,         &
-               rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av,             &
-               rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av
+        ONLY:  radiation, radiation_data_output_3d
 
 
@@ -494,60 +491,4 @@
              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_sw_cs_hr' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => rad_sw_cs_hr
-             ELSE
-                to_be_resorted => rad_sw_cs_hr_av
-             ENDIF
-
-          CASE ( 'rad_sw_hr' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => rad_sw_hr
-             ELSE
-                to_be_resorted => rad_sw_hr_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 ( 'rad_lw_cs_hr' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => rad_lw_cs_hr
-             ELSE
-                to_be_resorted => rad_lw_cs_hr_av
-             ENDIF
-
-          CASE ( 'rad_lw_hr' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => rad_lw_hr
-             ELSE
-                to_be_resorted => rad_lw_hr_av
-             ENDIF
-
           CASE ( 'rho' )
              IF ( av == 0 )  THEN
@@ -613,4 +554,23 @@
 
                 CALL lsm_data_output_3d( av, do3d(av,if), found, local_pf )
+                resorted = .TRUE.
+
+!
+!--             If no soil model variable was found, re-allocate local_pf
+                IF ( .NOT. found )  THEN
+                   nzb_do = nzb
+                   nzt_do = nz_do3d
+
+                   DEALLOCATE ( local_pf )
+                   ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )                 
+                ENDIF
+
+             ENDIF
+
+!
+!--          Radiation quantity
+             IF ( .NOT. found  .AND.  radiation )  THEN
+                CALL radiation_data_output_3d( av, do3d(av,if), found,         &
+                                               local_pf )
                 resorted = .TRUE.
              ENDIF
Index: palm/trunk/SOURCE/data_output_mask.f90
===================================================================
--- palm/trunk/SOURCE/data_output_mask.f90	(revision 1975)
+++ palm/trunk/SOURCE/data_output_mask.f90	(revision 1976)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Output of radiation quantities is now done directly in the respective module
 ! 
 ! Former revisions:
@@ -127,9 +127,5 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  rad_lw_in, rad_lw_in_av, rad_lw_out, rad_lw_out_av,             &
-               rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av,         &
-               rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av,             &
-               rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av
-
+        ONLY:  radiation, radiation_data_output_mask
 
     IMPLICIT NONE
@@ -425,60 +421,4 @@
              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 ( 'rad_lw_cs_hr' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => rad_lw_cs_hr
-             ELSE
-                to_be_resorted => rad_lw_cs_hr_av
-             ENDIF
-
-          CASE ( 'rad_lw_hr' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => rad_lw_hr
-             ELSE
-                to_be_resorted => rad_lw_hr_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_sw_cs_hr' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => rad_sw_cs_hr
-             ELSE
-                to_be_resorted => rad_sw_cs_hr_av
-             ENDIF
-
-          CASE ( 'rad_sw_hr' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => rad_sw_hr
-             ELSE
-                to_be_resorted => rad_sw_hr_av
-             ENDIF
-
           CASE ( 'rho' )
              IF ( av == 0 )  THEN
@@ -531,7 +471,19 @@
 
           CASE DEFAULT
+
+!
+!--          Radiation quantity
+             IF ( radiation )  THEN
+                CALL radiation_data_output_mask(av, domask(mid,av,if), found,  &
+                                                local_pf )
+             ENDIF
+
 !
 !--          User defined quantity
-             CALL user_data_output_mask(av, domask(mid,av,if), found, local_pf )
+             IF ( .NOT. found )  THEN
+                CALL user_data_output_mask(av, domask(mid,av,if), found,       &
+                                           local_pf )
+             ENDIF
+
              resorted = .TRUE.
 
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 1975)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 1976)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Removed some unneeded __rrtmg preprocessor directives
 ! 
 ! Former revisions:
@@ -832,7 +832,5 @@
                 ENDIF
 #endif
-
              ENDIF
-
 !
 !--          Subgridscale fluxes at the top surface
@@ -1447,5 +1445,4 @@
 
           IF ( radiation_scheme == 'rrtmg' )  THEN
-#if defined ( __rrtmg )
              hom(:,1,106,sr) = sums(:,106)            ! rad_lw_cs_hr
              hom(:,1,107,sr) = sums(:,107)            ! rad_lw_hr
@@ -1457,8 +1454,5 @@
              hom(:,1,112,sr) = sums(:,112)            ! rrtm_asdif
              hom(:,1,113,sr) = sums(:,113)            ! rrtm_asdir
-#endif
           ENDIF
-
-
        ENDIF
 
@@ -1643,5 +1637,4 @@
           ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_sw_out
 
-#if defined ( __rrtmg )
           IF ( radiation_scheme == 'rrtmg' )  THEN
              ts_value(dots_rad+5,sr) = hom(nzb,1,110,sr)          ! rrtm_aldif
@@ -1650,5 +1643,4 @@
              ts_value(dots_rad+8,sr) = hom(nzb,1,113,sr)          ! rrtm_asdir
           ENDIF
-#endif
 
        ENDIF
@@ -3619,5 +3611,4 @@
           ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_sw_out
 
-#if defined ( __rrtmg )
           IF ( radiation_scheme == 'rrtmg' )  THEN
              ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr)          ! rrtm_aldif
@@ -3626,5 +3617,4 @@
              ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr)          ! rrtm_asdir
           ENDIF
-#endif
 
        ENDIF
Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 1975)
+++ palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 1976)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Parts of the code have been reformatted. Use of radiation model output is
+! generalized and simplified. Added more output quantities due to modularization
 ! 
 ! Former revisions:
@@ -178,13 +179,7 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  force_radiation_call, radiation_scheme, rad_net, rad_sw_in,     &
-               rad_lw_out, rad_lw_out_change_0, sigma_sb,                      &
-               unscheduled_radiation_calls
+        ONLY:  force_radiation_call, rad_net, rad_sw_in, rad_lw_out,           &
+               rad_lw_out_change_0, unscheduled_radiation_calls
         
-#if defined ( __rrtmg )
-    USE radiation_model_mod,                                                   &
-        ONLY:  rrtm_idrv
-#endif
-
     USE statistics,                                                            &
         ONLY:  hom, statistic_regions
@@ -636,5 +631,4 @@
     END INTERFACE lsm_read_restart_data
 
-
     INTERFACE lsm_last_actions
        MODULE PROCEDURE lsm_last_actions
@@ -648,200 +642,206 @@
 !> Check data output for land surface model
 !------------------------------------------------------------------------------!
-    SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k )
+ SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k )
  
  
-       USE control_parameters,                                                 &
-           ONLY: data_output, message_string
-
-       IMPLICIT NONE
-
-       CHARACTER (LEN=*) ::  unit     !< 
-       CHARACTER (LEN=*) ::  var !<
-
-       INTEGER(iwp) :: i
-       INTEGER(iwp) :: ilen   
-       INTEGER(iwp) :: k
-
-       SELECT CASE ( TRIM( var ) )
-
-          CASE ( 'm_soil' )
-             IF (  .NOT.  land_surface )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                         'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
-             ENDIF
-             unit = 'm3/m3'
+    USE control_parameters,                                                 &
+        ONLY:  data_output, message_string
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  unit     !< 
+    CHARACTER (LEN=*) ::  var !<
+
+    INTEGER(iwp) :: i
+    INTEGER(iwp) :: ilen   
+    INTEGER(iwp) :: k
+
+    SELECT CASE ( TRIM( var ) )
+
+       CASE ( 'm_soil' )
+          IF (  .NOT.  land_surface )  THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                      'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+          ENDIF
+          unit = 'm3/m3'
+           
+       CASE ( 't_soil' )
+          IF (  .NOT.  land_surface )  THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                      'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+          ENDIF
+          unit = 'K'   
              
-          CASE ( 't_soil' )
-             IF (  .NOT.  land_surface )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                         'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
-             ENDIF
-             unit = 'K'   
+       CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'm_liq_eb*',&
+              'qsws_eb*', 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*',  &
+              'r_a*', 'r_s*', 'shf_eb*' )
+          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
+             message_string = 'illegal value for data_output: "' //         &
+                              TRIM( var ) // '" & only 2d-horizontal ' //   &
+                              'cross sections are allowed for this value'
+             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                              'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT.  land_surface )  THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                              'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT.  land_surface )  THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                              'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                              'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT.  land_surface )  THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                              'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT.  land_surface )  THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                              'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT.  land_surface )  THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                              'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )   &
+          THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                              'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT.  land_surface ) &
+          THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                              'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )   &
+          THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                              'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT.  land_surface ) &
+          THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                              'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT.  land_surface ) &
+          THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                              'res land_surface = .TRUE.'
+             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( TRIM( var ) == 'lai*'   )  unit = 'none' 
+          IF ( TRIM( var ) == 'c_liq*' )  unit = 'none'
+          IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
+          IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
+          IF ( TRIM( var ) == 'ghf_eb*')  unit = 'W/m2'
+          IF ( TRIM( var ) == 'm_liq_eb*'     )  unit = 'm'
+          IF ( TRIM( var ) == 'qsws_eb*'      ) unit = 'W/m2'
+          IF ( TRIM( var ) == 'qsws_liq_eb*'  ) unit = 'W/m2'
+          IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
+          IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
+          IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
+          IF ( TRIM( var ) == 'r_s*')     unit = 's/m' 
+          IF ( TRIM( var ) == 'shf_eb*')  unit = 'W/m2'
              
-          CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'm_liq_eb*',&
-                 'qsws_eb*', 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*',  &
-                 'r_a*', 'r_s*', 'shf_eb*' )
-             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
-                message_string = 'illegal value for data_output: "' //         &
-                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
-                                 'cross sections are allowed for this value'
-                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT.  land_surface )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT.  land_surface )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT.  land_surface )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT.  land_surface )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT.  land_surface )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )   &
-             THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT.  land_surface ) &
-             THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )   &
-             THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT.  land_surface ) &
-             THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
-             ENDIF
-             IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT.  land_surface ) &
-             THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE.'
-                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
-             ENDIF
-
-             IF ( TRIM( var ) == 'lai*'   )  unit = 'none' 
-             IF ( TRIM( var ) == 'c_liq*' )  unit = 'none'
-             IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
-             IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
-             IF ( TRIM( var ) == 'ghf_eb*')  unit = 'W/m2'
-             IF ( TRIM( var ) == 'm_liq_eb*'     )  unit = 'm'
-             IF ( TRIM( var ) == 'qsws_eb*'      ) unit = 'W/m2'
-             IF ( TRIM( var ) == 'qsws_liq_eb*'  ) unit = 'W/m2'
-             IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
-             IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
-             IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
-             IF ( TRIM( var ) == 'r_s*')     unit = 's/m' 
-             IF ( TRIM( var ) == 'shf_eb*')  unit = 'W/m2'
-             
-          CASE DEFAULT
-             unit = 'illegal'
-
-       END SELECT
-
-
-    END SUBROUTINE lsm_check_data_output
-
+       CASE DEFAULT
+          unit = 'illegal'
+
+    END SELECT
+
+
+ END SUBROUTINE lsm_check_data_output
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check data output of profiles for land surface model
+!------------------------------------------------------------------------------!
  SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit )
  
-       USE control_parameters,                                                 &
-           ONLY: data_output_pr, message_string
-
-       USE indices
-
-       USE profil_parameter
-
-       USE statistics
-
-       IMPLICIT NONE
+    USE control_parameters,                                                 &
+        ONLY:  data_output_pr, message_string
+
+    USE indices
+
+    USE profil_parameter
+
+    USE statistics
+
+    IMPLICIT NONE
    
-       CHARACTER (LEN=*) ::  unit      !< 
-       CHARACTER (LEN=*) ::  variable  !< 
-       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
+    CHARACTER (LEN=*) ::  unit      !< 
+    CHARACTER (LEN=*) ::  variable  !< 
+    CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
  
-       INTEGER(iwp) ::  user_pr_index !< 
-       INTEGER(iwp) ::  var_count     !< 
-
-       SELECT CASE ( TRIM( variable ) )
+    INTEGER(iwp) ::  user_pr_index !< 
+    INTEGER(iwp) ::  var_count     !< 
+
+    SELECT CASE ( TRIM( variable ) )
        
-          CASE ( 't_soil', '#t_soil' )
-             IF (  .NOT.  land_surface )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(var_count) ) // ' is' // &
-                                 'not implemented for land_surface = .FALSE.'
-                CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
-             ELSE
-                dopr_index(var_count) = 89
-                dopr_unit     = 'K'
-                hom(0:nzs-1,2,89,:)  = SPREAD( - zs, 2, statistic_regions+1 )
-                IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
-                   dopr_initial_index(var_count) = 90
-                   hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
-                   data_output_pr(var_count)     = data_output_pr(var_count)(2:)
-                ENDIF
-                unit = dopr_unit
+       CASE ( 't_soil', '#t_soil' )
+          IF (  .NOT.  land_surface )  THEN
+             message_string = 'data_output_pr = ' //                        &
+                              TRIM( data_output_pr(var_count) ) // ' is' // &
+                              'not implemented for land_surface = .FALSE.'
+             CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
+          ELSE
+             dopr_index(var_count) = 89
+             dopr_unit     = 'K'
+             hom(0:nzs-1,2,89,:)  = SPREAD( - zs, 2, statistic_regions+1 )
+             IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
+                dopr_initial_index(var_count) = 90
+                hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
+                data_output_pr(var_count)     = data_output_pr(var_count)(2:)
              ENDIF
-
-          CASE ( 'm_soil', '#m_soil' )
-             IF (  .NOT.  land_surface )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(var_count) ) // ' is' // &
-                                 ' not implemented for land_surface = .FALSE.'
-                CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
-             ELSE
-                dopr_index(var_count) = 91
-                dopr_unit     = 'm3/m3'
-                hom(0:nzs-1,2,91,:)  = SPREAD( - zs, 2, statistic_regions+1 )
-                IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
-                   dopr_initial_index(var_count) = 92
-                   hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
-                   data_output_pr(var_count)     = data_output_pr(var_count)(2:)
-                ENDIF
-                 unit = dopr_unit
+             unit = dopr_unit
+          ENDIF
+
+       CASE ( 'm_soil', '#m_soil' )
+          IF (  .NOT.  land_surface )  THEN
+             message_string = 'data_output_pr = ' //                        &
+                              TRIM( data_output_pr(var_count) ) // ' is' // &
+                              ' not implemented for land_surface = .FALSE.'
+             CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
+          ELSE
+             dopr_index(var_count) = 91
+             dopr_unit     = 'm3/m3'
+             hom(0:nzs-1,2,91,:)  = SPREAD( - zs, 2, statistic_regions+1 )
+             IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
+                dopr_initial_index(var_count) = 92
+                hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
+                data_output_pr(var_count)     = data_output_pr(var_count)(2:)
              ENDIF
-
-
-          CASE DEFAULT
-             unit = 'illegal'
-
-       END SELECT
-
-
-    END SUBROUTINE lsm_check_data_output_pr
+             unit = dopr_unit
+          ENDIF
+
+
+       CASE DEFAULT
+          unit = 'illegal'
+
+    END SELECT
+
+
+ END SUBROUTINE lsm_check_data_output_pr
  
  
@@ -851,190 +851,190 @@
 !> Check parameters routine for land surface model
 !------------------------------------------------------------------------------!
-    SUBROUTINE lsm_check_parameters
-
-       USE control_parameters,                                                 &
-           ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string,         &
-                 most_method, topography
+ SUBROUTINE lsm_check_parameters
+
+    USE control_parameters,                                                    &
+        ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string,           &
+               most_method, topography
                  
-       USE radiation_model_mod,                                                &
-           ONLY: radiation
+    USE radiation_model_mod,                                                   &
+        ONLY:  radiation
     
     
-       IMPLICIT NONE
+    IMPLICIT NONE
 
  
 !
-!--    Dirichlet boundary conditions are required as the surface fluxes are
-!--    calculated from the temperature/humidity gradients in the land surface
-!--    model
-       IF ( bc_pt_b == 'neumann'  .OR.  bc_q_b == 'neumann' )  THEN
-          message_string = 'lsm requires setting of'//                         &
-                           'bc_pt_b = "dirichlet" and '//                      &
-                           'bc_q_b  = "dirichlet"'
-          CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
+!-- Dirichlet boundary conditions are required as the surface fluxes are
+!-- calculated from the temperature/humidity gradients in the land surface
+!-- model
+    IF ( bc_pt_b == 'neumann'  .OR.  bc_q_b == 'neumann' )  THEN
+       message_string = 'lsm requires setting of'//                         &
+                        'bc_pt_b = "dirichlet" and '//                      &
+                        'bc_q_b  = "dirichlet"'
+       CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
+    ENDIF
+
+    IF (  .NOT.  constant_flux_layer )  THEN
+       message_string = 'lsm requires '//                                   &
+                        'constant_flux_layer = .T.'
+       CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
+    ENDIF
+
+    IF ( topography /= 'flat' )  THEN
+       message_string = 'lsm cannot be used ' //                            & 
+                        'in combination with  topography /= "flat"'
+       CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
+    ENDIF
+
+    IF ( ( veg_type == 14  .OR.  veg_type == 15 ) .AND.                     &
+           most_method == 'lookup' )  THEN
+        WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ',     &
+                                   'allowed in combination with ',          &
+                                   'most_method = ', most_method
+       CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
+    ENDIF
+
+    IF ( veg_type == 0 )  THEN
+       IF ( SUM( root_fraction ) /= 1.0_wp )  THEN
+          message_string = 'veg_type = 0 (user_defined)'//                  &
+                           'requires setting of root_fraction(0:3)'//       &
+                           '/= 9999999.9 and SUM(root_fraction) = 1'
+          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
        ENDIF
-
-       IF (  .NOT.  constant_flux_layer )  THEN
-          message_string = 'lsm requires '//                                   &
-                           'constant_flux_layer = .T.'
-          CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
+ 
+       IF ( min_canopy_resistance == 9999999.9_wp )  THEN
+          message_string = 'veg_type = 0 (user defined)'//                  &
+                           'requires setting of min_canopy_resistance'//    &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
        ENDIF
 
-       IF ( topography /= 'flat' )  THEN
-          message_string = 'lsm cannot be used ' //                            & 
-                           'in combination with  topography /= "flat"'
-          CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
+       IF ( leaf_area_index == 9999999.9_wp )  THEN
+          message_string = 'veg_type = 0 (user_defined)'//                  &
+                           'requires setting of leaf_area_index'//          &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
        ENDIF
 
-       IF ( ( veg_type == 14  .OR.  veg_type == 15 ) .AND.                     &
-              most_method == 'lookup' )  THEN
-           WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ',     &
-                                      'allowed in combination with ',          &
-                                      'most_method = ', most_method
-          CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
+       IF ( vegetation_coverage == 9999999.9_wp )  THEN
+          message_string = 'veg_type = 0 (user_defined)'//                  &
+                           'requires setting of vegetation_coverage'//      &
+                           '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
        ENDIF
 
-       IF ( veg_type == 0 )  THEN
-          IF ( SUM( root_fraction ) /= 1.0_wp )  THEN
-             message_string = 'veg_type = 0 (user_defined)'//                  &
-                              'requires setting of root_fraction(0:3)'//       &
-                              '/= 9999999.9 and SUM(root_fraction) = 1'
-             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-          ENDIF
- 
-          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
-             message_string = 'veg_type = 0 (user defined)'//                  &
-                              'requires setting of min_canopy_resistance'//    &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( leaf_area_index == 9999999.9_wp )  THEN
-             message_string = 'veg_type = 0 (user_defined)'//                  &
-                              'requires setting of leaf_area_index'//          &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( vegetation_coverage == 9999999.9_wp )  THEN
-             message_string = 'veg_type = 0 (user_defined)'//                  &
-                              'requires setting of vegetation_coverage'//      &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
-             message_string = 'veg_type = 0 (user_defined)'//                  &
-                              'requires setting of'//                          &
-                              'canopy_resistance_coefficient /= 9999999.9'
-             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
-             message_string = 'veg_type = 0 (user_defined)'//                  &
-                              'requires setting of lambda_surface_stable'//    &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
-             message_string = 'veg_type = 0 (user_defined)'//                  &
-                              'requires setting of lambda_surface_unstable'//  &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
-             message_string = 'veg_type = 0 (user_defined)'//                  &
-                              'requires setting of f_shortwave_incoming'//     &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( z0_eb == 9999999.9_wp )  THEN
-             message_string = 'veg_type = 0 (user_defined)'//                  &
-                              'requires setting of z0_eb'//                    &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( z0h_eb == 9999999.9_wp )  THEN
-             message_string = 'veg_type = 0 (user_defined)'//                  &
-                              'requires setting of z0h_eb'//                   &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
-          ENDIF
-
-
+       IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
+          message_string = 'veg_type = 0 (user_defined)'//                  &
+                           'requires setting of'//                          &
+                           'canopy_resistance_coefficient /= 9999999.9'
+          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
        ENDIF
 
-       IF ( soil_type == 0 )  THEN
-
-          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
-             message_string = 'soil_type = 0 (user_defined)'//                 &
-                              'requires setting of alpha_vangenuchten'//       &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( l_vangenuchten == 9999999.9_wp )  THEN
-             message_string = 'soil_type = 0 (user_defined)'//                 &
-                              'requires setting of l_vangenuchten'//           &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( n_vangenuchten == 9999999.9_wp )  THEN
-             message_string = 'soil_type = 0 (user_defined)'//                 &
-                              'requires setting of n_vangenuchten'//           &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
-             message_string = 'soil_type = 0 (user_defined)'//                 &
-                              'requires setting of hydraulic_conductivity'//   &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( saturation_moisture == 9999999.9_wp )  THEN
-             message_string = 'soil_type = 0 (user_defined)'//                 &
-                              'requires setting of saturation_moisture'//      &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( field_capacity == 9999999.9_wp )  THEN
-             message_string = 'soil_type = 0 (user_defined)'//                 &
-                              'requires setting of field_capacity'//           &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( wilting_point == 9999999.9_wp )  THEN
-             message_string = 'soil_type = 0 (user_defined)'//                 &
-                              'requires setting of wilting_point'//            &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          IF ( residual_moisture == 9999999.9_wp )  THEN
-             message_string = 'soil_type = 0 (user_defined)'//                 &
-                              'requires setting of residual_moisture'//        &
-                              '/= 9999999.9'
-             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
-          ENDIF
-
+       IF ( lambda_surface_stable == 9999999.9_wp )  THEN
+          message_string = 'veg_type = 0 (user_defined)'//                  &
+                           'requires setting of lambda_surface_stable'//    &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
        ENDIF
 
-       IF (  .NOT.  radiation )  THEN
-          message_string = 'lsm requires '//                                   &
-                           'radiation = .T.'
-          CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
+       IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
+          message_string = 'veg_type = 0 (user_defined)'//                  &
+                           'requires setting of lambda_surface_unstable'//  &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
        ENDIF
+
+       IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
+          message_string = 'veg_type = 0 (user_defined)'//                  &
+                           'requires setting of f_shortwave_incoming'//     &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( z0_eb == 9999999.9_wp )  THEN
+          message_string = 'veg_type = 0 (user_defined)'//                  &
+                           'requires setting of z0_eb'//                    &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( z0h_eb == 9999999.9_wp )  THEN
+          message_string = 'veg_type = 0 (user_defined)'//                  &
+                           'requires setting of z0h_eb'//                   &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+       ENDIF
+
+
+    ENDIF
+
+    IF ( soil_type == 0 )  THEN
+
+       IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
+          message_string = 'soil_type = 0 (user_defined)'//                 &
+                           'requires setting of alpha_vangenuchten'//       &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( l_vangenuchten == 9999999.9_wp )  THEN
+          message_string = 'soil_type = 0 (user_defined)'//                 &
+                           'requires setting of l_vangenuchten'//           &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( n_vangenuchten == 9999999.9_wp )  THEN
+          message_string = 'soil_type = 0 (user_defined)'//                 &
+                           'requires setting of n_vangenuchten'//           &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
+          message_string = 'soil_type = 0 (user_defined)'//                 &
+                           'requires setting of hydraulic_conductivity'//   &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( saturation_moisture == 9999999.9_wp )  THEN
+          message_string = 'soil_type = 0 (user_defined)'//                 &
+                           'requires setting of saturation_moisture'//      &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( field_capacity == 9999999.9_wp )  THEN
+          message_string = 'soil_type = 0 (user_defined)'//                 &
+                           'requires setting of field_capacity'//           &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( wilting_point == 9999999.9_wp )  THEN
+          message_string = 'soil_type = 0 (user_defined)'//                 &
+                           'requires setting of wilting_point'//            &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( residual_moisture == 9999999.9_wp )  THEN
+          message_string = 'soil_type = 0 (user_defined)'//                 &
+                           'requires setting of residual_moisture'//        &
+                           '/= 9999999.9'
+          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+       ENDIF
+
+    ENDIF
+
+    IF (  .NOT.  radiation )  THEN
+       message_string = 'lsm requires '//                                   &
+                        'radiation = .T.'
+       CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
+    ENDIF
        
        
-    END SUBROUTINE lsm_check_parameters
+ END SUBROUTINE lsm_check_parameters
  
 !------------------------------------------------------------------------------!
@@ -1043,504 +1043,463 @@
 !> Solver for the energy balance at the surface.
 !------------------------------------------------------------------------------!
-    SUBROUTINE lsm_energy_balance
-
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i         !< running index
-       INTEGER(iwp) ::  j         !< running index
-       INTEGER(iwp) ::  k, ks     !< running index
-
-       REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
-                   f1,          & !< resistance correction term 1
-                   f2,          & !< resistance correction term 2
-                   f3,          & !< resistance correction term 3
-                   m_min,       & !< minimum soil moisture
-                   e,           & !< water vapour pressure
-                   e_s,         & !< water vapour saturation pressure
-                   e_s_dt,      & !< derivate of e_s with respect to T
-                   tend,        & !< tendency
-                   dq_s_dt,     & !< derivate of q_s with respect to T
-                   coef_1,      & !< coef. for prognostic equation
-                   coef_2,      & !< coef. for prognostic equation
-                   f_qsws,      & !< factor for qsws_eb
-                   f_qsws_veg,  & !< factor for qsws_veg_eb
-                   f_qsws_soil, & !< factor for qsws_soil_eb
-                   f_qsws_liq,  & !< factor for qsws_liq_eb
-                   f_shf,       & !< factor for shf_eb
-                   lambda_surface, & !< Current value of lambda_surface
-                   m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
-                   pt1,         & !< potential temperature at first grid level
-                   qv1            !< specific humidity at first grid level
-
-!
-!--    Calculate the exner function for the current time step
-       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
-
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-             k = nzb_s_inner(j,i)
-
-!
-!--          Set lambda_surface according to stratification between skin layer and soil
-             IF (  .NOT.  pave_surface(j,i) )  THEN
-
-                c_surface_tmp = c_surface
-
-                IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i))  THEN
-                   lambda_surface = lambda_surface_s(j,i)
-                ELSE
-                   lambda_surface = lambda_surface_u(j,i)
-                ENDIF
+ SUBROUTINE lsm_energy_balance
+
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i         !< running index
+    INTEGER(iwp) ::  j         !< running index
+    INTEGER(iwp) ::  k, ks     !< running index
+
+    REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
+                f1,          & !< resistance correction term 1
+                f2,          & !< resistance correction term 2
+                f3,          & !< resistance correction term 3
+                m_min,       & !< minimum soil moisture
+                e,           & !< water vapour pressure
+                e_s,         & !< water vapour saturation pressure
+                e_s_dt,      & !< derivate of e_s with respect to T
+                tend,        & !< tendency
+                dq_s_dt,     & !< derivate of q_s with respect to T
+                coef_1,      & !< coef. for prognostic equation
+                coef_2,      & !< coef. for prognostic equation
+                f_qsws,      & !< factor for qsws_eb
+                f_qsws_veg,  & !< factor for qsws_veg_eb
+                f_qsws_soil, & !< factor for qsws_soil_eb
+                f_qsws_liq,  & !< factor for qsws_liq_eb
+                f_shf,       & !< factor for shf_eb
+                lambda_surface, & !< Current value of lambda_surface
+                m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
+                pt1,         & !< potential temperature at first grid level
+                qv1            !< specific humidity at first grid level
+
+!
+!-- Calculate the exner function for the current time step
+    exn = ( surface_pressure / 1000.0_wp )**0.286_wp
+
+    DO  i = nxlg, nxrg
+       DO  j = nysg, nyng
+          k = nzb_s_inner(j,i)
+
+!
+!--       Set lambda_surface according to stratification between skin layer and soil
+          IF (  .NOT.  pave_surface(j,i) )  THEN
+
+             c_surface_tmp = c_surface
+
+             IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i))  THEN
+                lambda_surface = lambda_surface_s(j,i)
              ELSE
-
-                c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
-                lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
-
+                lambda_surface = lambda_surface_u(j,i)
              ENDIF
-
-!
-!--          First step: calculate aerodyamic resistance. As pt, us, ts
-!--          are not available for the prognostic time step, data from the last
-!--          time step is used here. Note that this formulation is the
-!--          equivalent to the ECMWF formulation using drag coefficients
-             IF ( cloud_physics )  THEN
-                pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
-                qv1 = q(k+1,j,i) - ql(k+1,j,i)
-             ELSE
-                pt1 = pt(k+1,j,i)
-                qv1 = q(k+1,j,i)
-             ENDIF
-
-             r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
-
-!
-!--          Make sure that the resistance does not drop to zero for neutral 
-!--          stratification
-             IF ( ABS(r_a(j,i)) < 1.0_wp )  r_a(j,i) = 1.0_wp
-
-!
-!--          Second step: calculate canopy resistance r_canopy
-!--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
+          ELSE
+
+             c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
+             lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
+
+          ENDIF
+
+!
+!--       First step: calculate aerodyamic resistance. As pt, us, ts
+!--       are not available for the prognostic time step, data from the last
+!--       time step is used here. Note that this formulation is the
+!--       equivalent to the ECMWF formulation using drag coefficients
+          IF ( cloud_physics )  THEN
+             pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
+             qv1 = q(k+1,j,i) - ql(k+1,j,i)
+          ELSE
+             pt1 = pt(k+1,j,i)
+             qv1 = q(k+1,j,i)
+          ENDIF
+
+          r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
+
+!
+!--       Make sure that the resistance does not drop to zero for neutral 
+!--       stratification
+          IF ( ABS(r_a(j,i)) < 1.0_wp )  r_a(j,i) = 1.0_wp
+
+!
+!--       Second step: calculate canopy resistance r_canopy
+!--       f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
  
-!--          f1: correction for incoming shortwave radiation (stomata close at 
-!--          night)
-             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
-             ENDIF
-
-
-!
-!--          f2: correction for soil moisture availability to plants (the 
-!--          integrated soil moisture must thus be considered here)
-!--          f2 = 0 for very dry soils
-             m_total = 0.0_wp
-             DO  ks = nzb_soil, nzt_soil
-                 m_total = m_total + root_fr(ks,j,i)                           &
-                           * MAX(m_soil(ks,j,i),m_wilt(j,i))
-             ENDDO 
-
-             IF ( m_total > m_wilt(j,i)  .AND.  m_total < m_fc(j,i) )  THEN
-                f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
-             ELSEIF ( m_total >= m_fc(j,i) )  THEN
-                f2 = 1.0_wp
-             ELSE
-                f2 = 1.0E-20_wp
-             ENDIF
-
-!
-!--          Calculate water vapour pressure at saturation 
-             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
-                           - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
-
-!
-!--          f3: correction for vapour pressure deficit
-             IF ( g_d(j,i) /= 0.0_wp )  THEN
-!
-!--             Calculate vapour pressure
-                e  = qv1 * surface_pressure / 0.622_wp
-                f3 = EXP ( -g_d(j,i) * (e_s - e) )
-             ELSE
-                f3 = 1.0_wp
-             ENDIF
-
-!
-!--          Calculate canopy resistance. In case that c_veg is 0 (bare soils),
-!--          this calculation is obsolete, as r_canopy is not used below.
-!--          To do: check for very dry soil -> r_canopy goes to infinity
-             r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3      &
-                                             + 1.0E-20_wp)
-
-!
-!--          Third step: calculate bare soil resistance r_soil. The Clapp & 
-!--          Hornberger parametrization does not consider c_veg.
-             IF ( soil_type_2d(j,i) /= 7 )  THEN
-                m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *     &
-                        m_res(j,i)
-             ELSE
-                m_min = m_wilt(j,i)
-             ENDIF
-
-             f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
-             f2 = MAX(f2,1.0E-20_wp)
-             f2 = MIN(f2,1.0_wp)
-
-             r_soil(j,i) = r_soil_min(j,i) / f2
-
-!
-!--          Calculate the maximum possible liquid water amount on plants and
-!--          bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
-!--          assumed, while paved surfaces might hold up 1 mm of water. The 
-!--          liquid water fraction for paved surfaces is calculated after 
-!--          Noilhan & Planton (1989), while the ECMWF formulation is used for
-!--          vegetated surfaces and bare soils.
-             IF ( pave_surface(j,i) )  THEN
-                m_liq_eb_max = m_max_depth * 5.0_wp
-                c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 )
-             ELSE
-                m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i)           &
-                               + (1.0_wp - c_veg(j,i)) )
-                c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max )
-             ENDIF
-
-!
-!--          Calculate saturation specific humidity
-             q_s = 0.622_wp * e_s / surface_pressure
-
-!
-!--          In case of dewfall, set evapotranspiration to zero
-!--          All super-saturated water is then removed from the air
-             IF ( humidity  .AND.  q_s <= qv1 )  THEN
-                r_canopy(j,i) = 0.0_wp
-                r_soil(j,i)   = 0.0_wp
-             ENDIF
-
-!
-!--          Calculate coefficients for the total evapotranspiration 
-!--          In case of water surface, set vegetation and soil fluxes to zero.
-!--          For pavements, only evaporation of liquid water is possible.
-             IF ( water_surface(j,i) )  THEN
-                f_qsws_veg  = 0.0_wp
-                f_qsws_soil = 0.0_wp
-                f_qsws_liq  = rho_lv / r_a(j,i)
-             ELSEIF ( pave_surface (j,i) )  THEN
-                f_qsws_veg  = 0.0_wp
-                f_qsws_soil = 0.0_wp
-                f_qsws_liq  = rho_lv * c_liq(j,i) / r_a(j,i)
-             ELSE
-                f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/     &
-                              (r_a(j,i) + r_canopy(j,i))
-                f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +     &
-                                                          r_soil(j,i))
-                f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
-             ENDIF
-!
-!--          If soil moisture is below wilting point, plants do no longer
-!--          transpirate.
-!              IF ( m_soil(k,j,i) < m_wilt(j,i) )  THEN
-!                 f_qsws_veg = 0.0_wp
-!              ENDIF
-
-             f_shf  = rho_cp / r_a(j,i)
-             f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
-
-!
-!--          Calculate derivative of q_s for Taylor series expansion
-             e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -        &
-                              17.269_wp*(t_surface(j,i) - 273.16_wp)           &
-                              / (t_surface(j,i) - 35.86_wp)**2 )
-
-             dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
-
-!
-!--          Add LW up so that it can be removed in prognostic equation
-             rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
-
-!
-!--          Calculate new skin temperature
-             IF ( humidity )  THEN
-#if defined ( __rrtmg )
-!
-!--             Numerator of the prognostic equation
-                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
-                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
-                         + f_shf * pt1 + f_qsws * ( qv1 - q_s                  &
-                         + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
-                         * t_soil(nzb_soil,j,i)
-
-!
-!--             Denominator of the prognostic equation
-                coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt           &
-                         + lambda_surface + f_shf / exn
-#else
-
-!
-!--             Numerator of the prognostic equation
-                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
-                         * t_surface(j,i) ** 4                                 &
-                         + f_shf * pt1 + f_qsws * ( qv1                        &
-                         - q_s + dq_s_dt * t_surface(j,i) )                    &
-                         + lambda_surface * t_soil(nzb_soil,j,i)
-
-!
-!--             Denominator of the prognostic equation
-                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
-                         * dq_s_dt + lambda_surface + f_shf / exn
- 
-#endif
-             ELSE
-
-#if defined ( __rrtmg )
-!
-!--             Numerator of the prognostic equation
-                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
-                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
-                         + f_shf * pt1  + lambda_surface                       &
-                         * t_soil(nzb_soil,j,i)
-
-!
-!--             Denominator of the prognostic equation
-                coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
-#else
-
-!
-!--             Numerator of the prognostic equation
-                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
-                          * t_surface(j,i) ** 4 + f_shf * pt1                  &
-                          + lambda_surface * t_soil(nzb_soil,j,i)
-
-!
-!--             Denominator of the prognostic equation
-                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3               &
-                         + lambda_surface + f_shf / exn
-#endif
-             ENDIF
-
-             tend = 0.0_wp
-
-!
-!--          Implicit solution when the surface layer has no heat capacity,
-!--          otherwise use RK3 scheme.
-             t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *    &
-                                t_surface(j,i) ) / ( c_surface_tmp + coef_2    &
+!--       f1: correction for incoming shortwave radiation (stomata close at 
+!--       night)
+          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)) )
+
+
+
+!
+!--       f2: correction for soil moisture availability to plants (the 
+!--       integrated soil moisture must thus be considered here)
+!--       f2 = 0 for very dry soils
+          m_total = 0.0_wp
+          DO  ks = nzb_soil, nzt_soil
+              m_total = m_total + root_fr(ks,j,i)                              &
+                        * MAX(m_soil(ks,j,i),m_wilt(j,i))
+          ENDDO 
+
+          IF ( m_total > m_wilt(j,i)  .AND.  m_total < m_fc(j,i) )  THEN
+             f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
+          ELSEIF ( m_total >= m_fc(j,i) )  THEN
+             f2 = 1.0_wp
+          ELSE
+             f2 = 1.0E-20_wp
+          ENDIF
+
+!
+!--       Calculate water vapour pressure at saturation 
+          e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)        &
+                        - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
+
+!
+!--       f3: correction for vapour pressure deficit
+          IF ( g_d(j,i) /= 0.0_wp )  THEN
+!
+!--          Calculate vapour pressure
+             e  = qv1 * surface_pressure / 0.622_wp
+             f3 = EXP ( -g_d(j,i) * (e_s - e) )
+          ELSE
+             f3 = 1.0_wp
+          ENDIF
+
+!
+!--       Calculate canopy resistance. In case that c_veg is 0 (bare soils),
+!--       this calculation is obsolete, as r_canopy is not used below.
+!--       To do: check for very dry soil -> r_canopy goes to infinity
+          r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3         &
+                                          + 1.0E-20_wp)
+
+!
+!--       Third step: calculate bare soil resistance r_soil. The Clapp & 
+!--       Hornberger parametrization does not consider c_veg.
+          IF ( soil_type_2d(j,i) /= 7 )  THEN
+             m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *        &
+                     m_res(j,i)
+          ELSE
+             m_min = m_wilt(j,i)
+          ENDIF
+
+          f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
+          f2 = MAX(f2,1.0E-20_wp)
+          f2 = MIN(f2,1.0_wp)
+
+          r_soil(j,i) = r_soil_min(j,i) / f2
+
+!
+!--       Calculate the maximum possible liquid water amount on plants and
+!--       bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
+!--       assumed, while paved surfaces might hold up 1 mm of water. The 
+!--       liquid water fraction for paved surfaces is calculated after 
+!--       Noilhan & Planton (1989), while the ECMWF formulation is used for
+!--       vegetated surfaces and bare soils.
+          IF ( pave_surface(j,i) )  THEN
+             m_liq_eb_max = m_max_depth * 5.0_wp
+             c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 )
+          ELSE
+             m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i)              &
+                            + (1.0_wp - c_veg(j,i)) )
+             c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max )
+          ENDIF
+
+!
+!--       Calculate saturation specific humidity
+          q_s = 0.622_wp * e_s / surface_pressure
+
+!
+!--       In case of dewfall, set evapotranspiration to zero
+!--       All super-saturated water is then removed from the air
+          IF ( humidity  .AND.  q_s <= qv1 )  THEN
+             r_canopy(j,i) = 0.0_wp
+             r_soil(j,i)   = 0.0_wp
+          ENDIF
+
+!
+!--       Calculate coefficients for the total evapotranspiration 
+!--       In case of water surface, set vegetation and soil fluxes to zero.
+!--       For pavements, only evaporation of liquid water is possible.
+          IF ( water_surface(j,i) )  THEN
+             f_qsws_veg  = 0.0_wp
+             f_qsws_soil = 0.0_wp
+             f_qsws_liq  = rho_lv / r_a(j,i)
+          ELSEIF ( pave_surface (j,i) )  THEN
+             f_qsws_veg  = 0.0_wp
+             f_qsws_soil = 0.0_wp
+             f_qsws_liq  = rho_lv * c_liq(j,i) / r_a(j,i)
+          ELSE
+             f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/        &
+                           (r_a(j,i) + r_canopy(j,i))
+             f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +        &
+                                                             r_soil(j,i))
+             f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
+          ENDIF
+!
+!--       If soil moisture is below wilting point, plants do no longer
+!--       transpirate.
+!           IF ( m_soil(k,j,i) < m_wilt(j,i) )  THEN
+!              f_qsws_veg = 0.0_wp
+!           ENDIF
+
+          f_shf  = rho_cp / r_a(j,i)
+          f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
+
+!
+!--       Calculate derivative of q_s for Taylor series expansion
+          e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -           &
+                           17.269_wp*(t_surface(j,i) - 273.16_wp)              &
+                           / (t_surface(j,i) - 35.86_wp)**2 )
+
+          dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
+
+!
+!--       Add LW up so that it can be removed in prognostic equation
+          rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
+
+!
+!--       Calculate new skin temperature
+          IF ( humidity )  THEN
+
+!
+!--          Numerator of the prognostic equation
+             coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)                &
+                      * t_surface(j,i) - rad_lw_out(nzb,j,i)                   &
+                      + f_shf * pt1 + f_qsws * ( qv1 - q_s                     &
+                      + dq_s_dt * t_surface(j,i) ) + lambda_surface            &
+                      * t_soil(nzb_soil,j,i)
+
+!
+!--          Denominator of the prognostic equation
+             coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt              &
+                      + lambda_surface + f_shf / exn
+          ELSE
+
+!
+!--          Numerator of the prognostic equation
+             coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)                &
+                      * t_surface(j,i) - rad_lw_out(nzb,j,i)                   &
+                      + f_shf * pt1  + lambda_surface                          &
+                      * t_soil(nzb_soil,j,i)
+
+!
+!--          Denominator of the prognostic equation
+             coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
+
+          ENDIF
+
+          tend = 0.0_wp
+
+!
+!--       Implicit solution when the surface layer has no heat capacity,
+!--       otherwise use RK3 scheme.
+          t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *       &
+                             t_surface(j,i) ) / ( c_surface_tmp + coef_2       &
                                 * dt_3d * tsc(2) ) 
 
 !
-!--          Add RK3 term
-             IF ( c_surface_tmp /= 0.0_wp )  THEN
-
-                t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)           &
-                                   * tt_surface_m(j,i)
-
-!
-!--             Calculate true tendency
-                tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)     &
-                       * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
-!
-!--             Calculate t_surface tendencies for the next Runge-Kutta step
-                IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                   IF ( intermediate_timestep_count == 1 )  THEN
-                      tt_surface_m(j,i) = tend
-                   ELSEIF ( intermediate_timestep_count <                      &
-                            intermediate_timestep_count_max )  THEN
-                      tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp        &
-                                          * tt_surface_m(j,i)
-                   ENDIF
+!--       Add RK3 term
+          IF ( c_surface_tmp /= 0.0_wp )  THEN
+
+             t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)              &
+                                * tt_surface_m(j,i)
+
+!
+!--          Calculate true tendency
+             tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)        &
+                    * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
+!
+!--          Calculate t_surface tendencies for the next Runge-Kutta step
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( intermediate_timestep_count == 1 )  THEN
+                   tt_surface_m(j,i) = tend
+                ELSEIF ( intermediate_timestep_count <                         &
+                         intermediate_timestep_count_max )  THEN
+                   tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp           &
+                                       * tt_surface_m(j,i)
                 ENDIF
              ENDIF
-
-!
-!--          In case of fast changes in the skin temperature, it is possible to
-!--          update the radiative fluxes independently from the prescribed
-!--          radiation call frequency. This effectively prevents oscillations,
-!--          especially when setting skip_time_do_radiation /= 0. The threshold
-!--          value of 0.2 used here is just a first guess. This method should be
-!--          revised in the future as tests have shown that the threshold is 
-!--          often reached, when no oscillations would occur (causes immense
-!--          computing time for the radiation code).
-             IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp  .AND.     &
-                  unscheduled_radiation_calls )  THEN
-                force_radiation_call_l = .TRUE.
+          ENDIF
+
+!
+!--       In case of fast changes in the skin temperature, it is possible to
+!--       update the radiative fluxes independently from the prescribed
+!--       radiation call frequency. This effectively prevents oscillations,
+!--       especially when setting skip_time_do_radiation /= 0. The threshold
+!--       value of 0.2 used here is just a first guess. This method should be
+!--       revised in the future as tests have shown that the threshold is 
+!--       often reached, when no oscillations would occur (causes immense
+!--       computing time for the radiation code).
+          IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp  .AND.        &
+               unscheduled_radiation_calls )  THEN
+             force_radiation_call_l = .TRUE.
+          ENDIF
+
+          pt(k,j,i) = t_surface_p(j,i) / exn
+
+!
+!--       Calculate fluxes
+          rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)         &
+                             * t_surface(j,i) - rad_lw_out(nzb,j,i)            &
+                             - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
+
+          rad_net(j,i) = rad_net_l(j,i)
+          rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) + rad_lw_out_change_0(j,i) &
+                                * ( t_surface_p(j,i) - t_surface(j,i) )
+
+          ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)                  &
+                           - t_soil(nzb_soil,j,i))
+
+          shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
+
+          shf(j,i) = shf_eb(j,i) / rho_cp
+
+          IF ( humidity )  THEN
+             qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt               &
+                             * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
+
+             qsws(j,i) = qsws_eb(j,i) / rho_lv
+
+             qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                   &
+                                 + dq_s_dt * t_surface(j,i) - dq_s_dt          &
+                                 * t_surface_p(j,i) )
+
+             qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                   &
+                                 + dq_s_dt * t_surface(j,i) - dq_s_dt          &
+                                 * t_surface_p(j,i) )
+
+             qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv1 - q_s                   &
+                                 + dq_s_dt * t_surface(j,i) - dq_s_dt          &
+                                 * t_surface_p(j,i) )
+          ENDIF
+
+!
+!--       Calculate the true surface resistance
+          IF ( qsws_eb(j,i) == 0.0_wp )  THEN
+             r_s(j,i) = 1.0E10_wp
+          ELSE
+             r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                       &
+                        * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )        &
+                        / qsws_eb(j,i) - r_a(j,i)
+          ENDIF
+
+!
+!--       Calculate change in liquid water reservoir due to dew fall or 
+!--       evaporation of liquid water
+          IF ( humidity )  THEN
+!
+!--          If precipitation is activated, add rain water to qsws_liq_eb
+!--          and qsws_soil_eb according the the vegetation coverage.
+!--          precipitation_rate is given in mm.
+             IF ( precipitation )  THEN
+
+!
+!--             Add precipitation to liquid water reservoir, if possible.
+!--             Otherwise, add the water to soil. In case of
+!--             pavements, the exceeding water amount is implicitely removed 
+!--             as runoff as qsws_soil_eb is then not used in the soil model
+                IF ( m_liq_eb(j,i) /= m_liq_eb_max )  THEN
+                   qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                         &
+                                    + c_veg(j,i) * prr(k,j,i) * hyrho(k)       &
+                                    * 0.001_wp * rho_l * l_v
+                ELSE
+                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                       &
+                                     + c_veg(j,i) * prr(k,j,i) * hyrho(k)      &
+                                     * 0.001_wp * rho_l * l_v
+                ENDIF
+
+!--             Add precipitation to bare soil according to the bare soil
+!--             coverage.
+                qsws_soil_eb(j,i) = qsws_soil_eb(j,i) + (1.0_wp                &
+                                    - c_veg(j,i)) * prr(k,j,i) * hyrho(k)      &
+                                    * 0.001_wp * rho_l * l_v
              ENDIF
 
-             pt(k,j,i) = t_surface_p(j,i) / exn
-
-!
-!--          Calculate fluxes
-#if defined ( __rrtmg )
-             rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)      &
-                                * t_surface(j,i) - rad_lw_out(nzb,j,i)         &
-                                - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
-
-             IF ( rrtm_idrv == 1 )  THEN
-                rad_net(j,i) = rad_net_l(j,i)
-                rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i)                      &
-                                      + rad_lw_out_change_0(j,i)               &
-                                      * ( t_surface_p(j,i) - t_surface(j,i) )
+!
+!--          If the air is saturated, check the reservoir water level
+             IF ( qsws_eb(j,i) < 0.0_wp )  THEN
+
+!
+!--             Check if reservoir is full (avoid values > m_liq_eb_max)
+!--             In that case, qsws_liq_eb goes to qsws_soil_eb. In this 
+!--             case qsws_veg_eb is zero anyway (because c_liq = 1),       
+!--             so that tend is zero and no further check is needed
+                IF ( m_liq_eb(j,i) == m_liq_eb_max )  THEN
+                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                       &
+                                        + qsws_liq_eb(j,i)
+
+                   qsws_liq_eb(j,i)  = 0.0_wp
+                ENDIF
+
+!
+!--             In case qsws_veg_eb becomes negative (unphysical behavior), 
+!--             let the water enter the liquid water reservoir as dew on the
+!--             plant
+                IF ( qsws_veg_eb(j,i) < 0.0_wp )  THEN
+                   qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
+                   qsws_veg_eb(j,i) = 0.0_wp
+                ENDIF
+             ENDIF                    
+  
+             tend = - qsws_liq_eb(j,i) * drho_l_lv
+             m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend         &
+                                                + tsc(3) * tm_liq_eb_m(j,i) )
+
+!
+!--          Check if reservoir is overfull -> reduce to maximum
+!--          (conservation of water is violated here)
+             m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
+
+!
+!--          Check if reservoir is empty (avoid values < 0.0)
+!--          (conservation of water is violated here)
+             m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
+
+
+!
+!--          Calculate m_liq_eb tendencies for the next Runge-Kutta step
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( intermediate_timestep_count == 1 )  THEN
+                   tm_liq_eb_m(j,i) = tend
+                ELSEIF ( intermediate_timestep_count <                         &
+                         intermediate_timestep_count_max )  THEN
+                   tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp            &
+                                    * tm_liq_eb_m(j,i)
+                ENDIF
              ENDIF
+
+          ENDIF
+
+       ENDDO
+    ENDDO
+
+!
+!-- Make a logical OR for all processes. Force radiation call if at
+!-- least one processor reached the threshold change in skin temperature
+    IF ( unscheduled_radiation_calls  .AND.  intermediate_timestep_count       &
+         == intermediate_timestep_count_max-1 )  THEN
+#if defined( __parallel )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
+                           1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
 #else
-             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)
+       force_radiation_call = force_radiation_call_l
 #endif
-
-             ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
-                              - t_soil(nzb_soil,j,i))
-
-             shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
-
-             shf(j,i) = shf_eb(j,i) / rho_cp
-
-             IF ( humidity )  THEN
-                qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt            &
-                                * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
-
-                qsws(j,i) = qsws_eb(j,i) / rho_lv
-
-                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                &
-                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
-                                    * t_surface_p(j,i) )
-
-                qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                &
-                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
-                                    * t_surface_p(j,i) )
-
-                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv1 - q_s                &
-                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
-                                    * t_surface_p(j,i) )
-             ENDIF
-
-!
-!--          Calculate the true surface resistance
-             IF ( qsws_eb(j,i) == 0.0_wp )  THEN
-                r_s(j,i) = 1.0E10_wp
-             ELSE
-                r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                    &
-                           * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
-                           / qsws_eb(j,i) - r_a(j,i)
-             ENDIF
-
-!
-!--          Calculate change in liquid water reservoir due to dew fall or 
-!--          evaporation of liquid water
-             IF ( humidity )  THEN
-!
-!--             If precipitation is activated, add rain water to qsws_liq_eb
-!--             and qsws_soil_eb according the the vegetation coverage.
-!--             precipitation_rate is given in mm.
-                IF ( precipitation )  THEN
-
-!
-!--                Add precipitation to liquid water reservoir, if possible.
-!--                Otherwise, add the water to soil. In case of
-!--                pavements, the exceeding water amount is implicitely removed 
-!--                as runoff as qsws_soil_eb is then not used in the soil model
-                   IF ( m_liq_eb(j,i) /= m_liq_eb_max )  THEN
-                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                      &
-                                       + c_veg(j,i) * prr(k,j,i) * hyrho(k)    &
-                                       * 0.001_wp * rho_l * l_v
-                   ELSE
-                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
-                                        + c_veg(j,i) * prr(k,j,i) * hyrho(k)   &
-                                        * 0.001_wp * rho_l * l_v
-                   ENDIF
-
-!--                Add precipitation to bare soil according to the bare soil
-!--                coverage.
-                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i) + (1.0_wp             &
-                                       - c_veg(j,i)) * prr(k,j,i) * hyrho(k)   &
-                                       * 0.001_wp * rho_l * l_v
-                ENDIF
-
-!
-!--             If the air is saturated, check the reservoir water level
-                IF ( qsws_eb(j,i) < 0.0_wp )  THEN
-
-!
-!--                Check if reservoir is full (avoid values > m_liq_eb_max)
-!--                In that case, qsws_liq_eb goes to qsws_soil_eb. In this 
-!--                case qsws_veg_eb is zero anyway (because c_liq = 1),       
-!--                so that tend is zero and no further check is needed
-                   IF ( m_liq_eb(j,i) == m_liq_eb_max )  THEN
-                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
-                                           + qsws_liq_eb(j,i)
-
-                      qsws_liq_eb(j,i)  = 0.0_wp
-                   ENDIF
-
-!
-!--                In case qsws_veg_eb becomes negative (unphysical behavior), 
-!--                let the water enter the liquid water reservoir as dew on the
-!--                plant
-                   IF ( qsws_veg_eb(j,i) < 0.0_wp )  THEN
-                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
-                      qsws_veg_eb(j,i) = 0.0_wp
-                   ENDIF
-                ENDIF                    
-  
-                tend = - qsws_liq_eb(j,i) * drho_l_lv
-
-                m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend      &
-                                                   + tsc(3) * tm_liq_eb_m(j,i) )
-
-!
-!--             Check if reservoir is overfull -> reduce to maximum
-!--             (conservation of water is violated here)
-                m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
-
-!
-!--             Check if reservoir is empty (avoid values < 0.0)
-!--             (conservation of water is violated here)
-                m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
-
-
-!
-!--             Calculate m_liq_eb tendencies for the next Runge-Kutta step
-                IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                   IF ( intermediate_timestep_count == 1 )  THEN
-                      tm_liq_eb_m(j,i) = tend
-                   ELSEIF ( intermediate_timestep_count <                      &
-                            intermediate_timestep_count_max )  THEN
-                      tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp         &
-                                      * tm_liq_eb_m(j,i)
-                   ENDIF
-                ENDIF
-
-             ENDIF
-
-          ENDDO
-       ENDDO
-
-!
-!--    Make a logical OR for all processes. Force radiation call if at
-!--    least one processor reached the threshold change in skin temperature
-       IF ( unscheduled_radiation_calls  .AND.  intermediate_timestep_count    &
-            == intermediate_timestep_count_max-1 )  THEN
-#if defined( __parallel )
-          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
-                              1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
-#else
-          force_radiation_call = force_radiation_call_l
-#endif
-          force_radiation_call_l = .FALSE.
-       ENDIF
-
-!
-!--    Calculate surface specific humidity
-       IF ( humidity )  THEN
-          CALL calc_q_surface
-       ENDIF
-
-!
-!--    Calculate new roughness lengths (for water surfaces only)
-       CALL calc_z0_water_surface
-
-
-    END SUBROUTINE lsm_energy_balance
+       force_radiation_call_l = .FALSE.
+    ENDIF
+
+!
+!-- Calculate surface specific humidity
+    IF ( humidity )  THEN
+       CALL calc_q_surface
+    ENDIF
+
+!
+!-- Calculate new roughness lengths (for water surfaces only)
+    CALL calc_z0_water_surface
+
+
+ END SUBROUTINE lsm_energy_balance
+
 
 !------------------------------------------------------------------------------!
@@ -2451,5 +2410,5 @@
 ! Description:
 ! ------------
-!> Soubroutine for averaging 3D data
+!> Subroutine for averaging 3D data
 !------------------------------------------------------------------------------!
 SUBROUTINE lsm_3d_data_averaging( mode, variable )
@@ -2805,37 +2764,37 @@
 ! Description:
 ! ------------
-!> Soubroutine defines appropriate grid for netcdf variables.
+!> Subroutine defining appropriate grid for netcdf variables.
 !> It is called out from subroutine netcdf.
 !------------------------------------------------------------------------------!
-    SUBROUTINE lsm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
+ SUBROUTINE lsm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
     
-        IMPLICIT NONE
-
-        CHARACTER (LEN=*), INTENT(IN)  ::  var         !< 
-        LOGICAL, INTENT(OUT)           ::  found       !< 
-        CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !< 
-        CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !< 
-        CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !< 
-
-
-
-!
-!--     Check for the grid
-        SELECT CASE ( TRIM( var ) )
-
-           CASE ( 'm_soil', 't_soil' )
-              found  = .TRUE.
-              grid_x = 'x'
-              grid_y = 'y'
-              grid_z = 'zs'
-
-           CASE DEFAULT
-              found  = .FALSE.
-              grid_x = 'none'
-              grid_y = 'none'
-              grid_z = 'none'
-        END SELECT
-
-    END SUBROUTINE lsm_define_netcdf_grid
+     IMPLICIT NONE
+
+     CHARACTER (LEN=*), INTENT(IN)  ::  var         !< 
+     LOGICAL, INTENT(OUT)           ::  found       !< 
+     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !< 
+     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !< 
+     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !< 
+
+     found  = .TRUE.
+
+!
+!--  Check for the grid
+     SELECT CASE ( TRIM( var ) )
+
+        CASE ( 'm_soil', 't_soil', 'm_soil_xy', 't_soil_xy', 'm_soil_xz',      &
+               't_soil_xz', 'm_soil_yz', 't_soil_yz' )
+           grid_x = 'x'
+           grid_y = 'y'
+           grid_z = 'zs'
+
+        CASE DEFAULT
+           found  = .FALSE.
+           grid_x = 'none'
+           grid_y = 'none'
+           grid_z = 'none'
+     END SELECT
+
+ END SUBROUTINE lsm_define_netcdf_grid
 
 !------------------------------------------------------------------------------!
@@ -2843,15 +2802,13 @@
 ! Description:
 ! ------------
-!> Soubroutine defines 3D output variables
+!> Subroutine defining 3D output variables
 !------------------------------------------------------------------------------!
  SUBROUTINE lsm_data_output_2d( av, variable, found, grid, mode, local_pf,     &
                                 two_d, nzb_do, nzt_do )
  
-
     USE indices
 
     USE kinds
 
-    USE user
 
     IMPLICIT NONE
@@ -3174,5 +3131,5 @@
 ! Description:
 ! ------------
-!> Soubroutine defines 3D output variables
+!> Subroutine defining 3D output variables
 !------------------------------------------------------------------------------!
  SUBROUTINE lsm_data_output_3d( av, variable, found, local_pf )
Index: palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 1975)
+++ palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 1976)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Removed remaining 2D land surface quantities. Definition of radiation 
+! quantities is now done directly in the respective module
 ! 
 ! Former revisions:
@@ -148,4 +149,6 @@
 !> @todo calculation of output time levels for parallel NetCDF still does not
 !>       cover every exception (change of dt_do, end_time in restart)
+!> @todo timeseries and profile output still needs to be rewritten to allow
+!>       modularization
 !------------------------------------------------------------------------------!
  MODULE netcdf_interface
@@ -414,4 +417,7 @@
     USE profil_parameter,                                                      &
         ONLY:  crmax, cross_profiles, dopr_index, profile_columns, profile_rows
+
+    USE radiation_model_mod,                                                   &
+        ONLY: radiation, radiation_define_netcdf_grid
 
     USE spectra_mod,                                                           &
@@ -735,6 +741,5 @@
                 CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
                        'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',        &
-                       'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr',            &
-                       'rad_sw_hr', 'rho', 's', 'sa', 'vpt'  )
+                       'rho', 's', 'sa', 'vpt'  )
 
                    grid_x = 'x'
@@ -757,6 +762,5 @@
 !
 !--             w grid
-                CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',   &
-                       'w' )
+                CASE ( 'w' )
 
                    grid_x = 'x'
@@ -772,4 +776,12 @@
                       CALL lsm_define_netcdf_grid( domask(mid,av,i), found,    &
                                                    grid_x, grid_y, grid_z )
+                   ENDIF
+
+!
+!--                Check for radiation quantities
+                   IF ( .NOT. found  .AND.  radiation )  THEN
+                      CALL radiation_define_netcdf_grid( domask(mid,av,i),     &
+                                                         found, grid_x, grid_y,&
+                                                         grid_z )
                    ENDIF
 
@@ -1239,6 +1251,5 @@
                 CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
                        'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho', &
-                       's', 'sa', 'vpt' , 'rad_lw_cs_hr', 'rad_lw_hr',         &
-                       'rad_sw_cs_hr', 'rad_sw_hr' )
+                       's', 'sa', 'vpt' )
 
                    grid_x = 'x'
@@ -1261,6 +1272,5 @@
 !
 !--             w grid
-                CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',   &
-                       'w' )
+                CASE ( 'w' )
 
                    grid_x = 'x'
@@ -1275,5 +1285,13 @@
                    IF ( land_surface )  THEN
                       CALL lsm_define_netcdf_grid( do3d(av,i), found, grid_x,  &
-                                                    grid_y, grid_z )
+                                                   grid_y, grid_z )
+                   ENDIF
+
+!
+!--                Check for radiation quantities
+                   IF ( .NOT. found  .AND.  radiation )  THEN
+                      CALL radiation_define_netcdf_grid( do3d(av,i), found,    &
+                                                         grid_x, grid_y,       &
+                                                         grid_z )
                    ENDIF
                    
@@ -1837,7 +1855,6 @@
                              'pr_xy', 'prr_xy', 'pt_xy', 'q_xy', 'qc_xy',      &
                              'ql_xy', 'ql_c_xy', 'ql_v_xy', 'ql_vp_xy',        &
-                             'qr_xy', 'qv_xy', 'rad_lw_cs_hr_xy',              &
-                             'rad_lw_hr_xy', 'rad_sw_cs_hr_xy', 'rad_sw_hr_xy',& 
-                             'rho_xy', 's_xy', 'sa_xy', 'vpt_xy' )
+                             'qr_xy', 'qv_xy', 'rho_xy', 's_xy', 'sa_xy',      &
+                             'vpt_xy' )
 
                          grid_x = 'x'
@@ -1860,22 +1877,34 @@
 !
 !--                   w grid
-                      CASE ( 'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy',  &
-                             'rad_sw_out_xy' , 'w_xy' )
+                      CASE ( 'w_xy' )
 
                          grid_x = 'x'
                          grid_y = 'y'
                          grid_z = 'zw'
-!
-!--                   soil grid
-                      CASE ( 'm_soil_xy', 't_soil_xy' )
-                         grid_x = 'x'
-                         grid_y = 'y'
-                         grid_z = 'zs'
+
 
                       CASE DEFAULT
 !
+!--                      Check for land surface quantities
+                         IF ( land_surface )  THEN
+                            CALL lsm_define_netcdf_grid( do2d(av,i), found,    &
+                                                   grid_x, grid_y, grid_z )
+                         ENDIF
+
+!
+!--                      Check for radiation quantities
+                         IF ( .NOT. found  .AND.  radiation )  THEN
+                            CALL radiation_define_netcdf_grid( do2d(av,i),     &
+                                                         found, grid_x, grid_y,&
+                                                         grid_z )
+                         ENDIF
+
+!
 !--                      Check for user-defined quantities
-                         CALL user_define_netcdf_grid( do2d(av,i), found, &
-                                                       grid_x, grid_y, grid_z )
+                         IF ( .NOT. found )  THEN
+                            CALL user_define_netcdf_grid( do2d(av,i), found,   &
+                                                          grid_x, grid_y,      &
+                                                          grid_z )
+                         ENDIF
 
                          IF ( .NOT. found )  THEN
@@ -2491,7 +2520,5 @@
                           'prr_xz', 'pt_xz', 'q_xz', 'qc_xz', 'ql_xz',         &
                           'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz', 'qv_xz',  &
-                          'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',                   &
-                          'rad_sw_cs_hr_xz', 'rad_sw_hr_xz''rho_xz', 's_xz',   &
-                          'sa_xz', 'vpt_xz' )
+                          'rho_xz', 's_xz', 'sa_xz', 'vpt_xz' )
 
                       grid_x = 'x'
@@ -2514,6 +2541,5 @@
 !
 !--                w grid
-                   CASE ( 'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz',     &
-                          'rad_sw_out_xz', 'w_xz' )
+                   CASE ( 'w_xz' )
 
                       grid_x = 'x'
@@ -2521,17 +2547,28 @@
                       grid_z = 'zw'
 
-!
-!--                soil grid
-                   CASE ( 'm_soil_xz', 't_soil_xz' )
-
-                      grid_x = 'x'
-                      grid_y = 'y'
-                      grid_z = 'zs'
-
                    CASE DEFAULT
+
+!
+!--                   Check for land surface quantities
+                      IF ( land_surface )  THEN
+                         CALL lsm_define_netcdf_grid( do2d(av,i), found,       &
+                                                      grid_x, grid_y, grid_z )
+                      ENDIF
+
+!
+!--                   Check for radiation quantities
+                      IF ( .NOT. found  .AND.  radiation )  THEN
+                         CALL radiation_define_netcdf_grid( do2d(av,i), found, &
+                                                            grid_x, grid_y,    &
+                                                            grid_z )
+                      ENDIF
+
 !
 !--                   Check for user-defined quantities
-                      CALL user_define_netcdf_grid( do2d(av,i), found, &
-                                                    grid_x, grid_y, grid_z )
+                      IF ( .NOT. found )  THEN
+                         CALL user_define_netcdf_grid( do2d(av,i), found,      &
+                                                       grid_x, grid_y, grid_z )
+                      ENDIF
+
                       IF ( .NOT. found )  THEN
                          WRITE ( message_string, * ) 'no grid defined for', &
@@ -3137,7 +3174,5 @@
                           'prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz',         &
                           'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz',  &
-                          'rad_lw_cs_hr_yz', 'rad_lw_hr_yz',                   &
-                          'rad_sw_cs_hr_yz', 'rad_sw_hr_yz''rho_yz', 's_yz',   &
-                          'sa_yz', 'vpt_yz' )
+                          'rho_yz', 's_yz', 'sa_yz', 'vpt_yz' )
 
                       grid_x = 'x'
@@ -3160,23 +3195,33 @@
 !
 !--                w grid
-                   CASE ( 'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz',     &
-                          'rad_sw_out_yz', 'w_yz' )
+                   CASE ( 'w_yz' )
 
                       grid_x = 'x'
                       grid_y = 'y'
                       grid_z = 'zw'
-!
-!--                soil grid
-                   CASE ( 'm_soil_yz', 't_soil_yz' )
-
-                      grid_x = 'x'
-                      grid_y = 'y'
-                      grid_z = 'zs'
+
 
                    CASE DEFAULT
 !
+!--                   Check for land surface quantities
+                      IF ( land_surface )  THEN
+                         CALL lsm_define_netcdf_grid( do2d(av,i), found,       &
+                                                      grid_x, grid_y, grid_z )
+                      ENDIF
+
+!
+!--                   Check for radiation quantities
+                      IF ( .NOT. found  .AND.  radiation )  THEN
+                         CALL radiation_define_netcdf_grid( do2d(av,i), found, &
+                                                            grid_x, grid_y,    &
+                                                            grid_z )
+                      ENDIF
+
+!
 !--                   Check for user-defined quantities
-                      CALL user_define_netcdf_grid( do2d(av,i), found, &
-                                                    grid_x, grid_y, grid_z )
+                      IF ( .NOT. found )  THEN
+                         CALL user_define_netcdf_grid( do2d(av,i), found,      &
+                                                       grid_x, grid_y, grid_z )
+                      ENDIF
 
                       IF ( .NOT. found )  THEN
@@ -4313,5 +4358,5 @@
 !--                Check for user-defined quantities (found, grid_x and grid_y
 !--                are dummies)
-                   CALL user_define_netcdf_grid( data_output_sp(i), found,  &
+                   CALL user_define_netcdf_grid( data_output_sp(i), found,     &
                                                  grid_x, grid_y, grid_z )
 
Index: palm/trunk/SOURCE/palm.f90
===================================================================
--- palm/trunk/SOURCE/palm.f90	(revision 1975)
+++ palm/trunk/SOURCE/palm.f90	(revision 1976)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added call to radiation_last_actions for binary output of land surface model 
+! data
 ! 
 ! Former revisions:
@@ -134,4 +135,6 @@
 !> machines are controlled via cpp-directives.
 !> Model runs are only feasible using the ksh-script mrun.
+!>
+!> @todo create routine last_actions instead of calling lsm_last_actions etc.
 !------------------------------------------------------------------------------!
  PROGRAM palm
@@ -182,4 +185,7 @@
         ONLY:  cpl_id, nested_run, pmci_child_initialize, pmci_init,           &
                pmci_modelconfiguration, pmci_parent_initialize
+
+    USE radiation_model_mod,                                                   &
+        ONLY:  radiation, radiation_last_actions
 
     USE statistics,                                                            &
@@ -376,5 +382,4 @@
     ENDIF
 
-
 !
 !-- Output of program header
@@ -394,4 +399,5 @@
        CALL data_output_2d( 'yz', 0 )
     ENDIF
+
     IF ( do3d_at_begin )  THEN
        CALL data_output_3d( 0 )
@@ -458,4 +464,7 @@
              CALL lsm_last_actions
           ENDIF
+          IF ( radiation )  THEN
+             CALL radiation_last_actions
+          ENDIF
           CALL user_last_actions
           IF ( write_binary(1:4) == 'true' )  CALL close_file( 14 )
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 1975)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 1976)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Simplied calls to radiation model
 ! 
 ! Former revisions:
@@ -305,5 +305,5 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  radiation, radiation_scheme, radiation_tendency,                &
+        ONLY:  radiation, radiation_tendency,                                  &
                skip_time_do_radiation
 
@@ -638,5 +638,5 @@
 !
 !--          If required, add tendency due to radiative heating/cooling
-             IF ( radiation .AND. radiation_scheme == 'rrtmg' .AND.            &
+             IF ( radiation  .AND.                                             &
                   simulated_time > skip_time_do_radiation )  THEN
                 CALL radiation_tendency ( i, j, tend )
@@ -1371,5 +1371,5 @@
 !
 !--    If required, add tendency due to radiative heating/cooling
-       IF ( radiation .AND. radiation_scheme == 'rrtmg' .AND.                  &
+       IF ( radiation  .AND.                                                   &
             simulated_time > skip_time_do_radiation )  THEN
             CALL radiation_tendency ( tend )
@@ -2293,5 +2293,5 @@
        ENDIF
 
-       IF ( radiation .AND. radiation_scheme == 'rrtmg' .AND.                  &
+       IF ( radiation .AND.                                                    &
             simulated_time > skip_time_do_radiation )  THEN
             CALL radiation_tendency ( tend )
Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 1975)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 1976)
@@ -19,5 +19,8 @@
 ! Current revisions:
 ! -----------------
-! 
+! Output of 2D/3D/masked data is now directly done within this module. The
+! radiation schemes have been simplified for better usability so that
+! rad_lw_in, rad_lw_out, rad_sw_in, and rad_sw_out are available independent of
+! the radiation code used.
 ! 
 ! Former revisions:
@@ -130,5 +133,4 @@
 
 #if defined ( __rrtmg )
-
     USE parrrsw,                                                               &
         ONLY:  naerec, nbndsw
@@ -276,5 +278,6 @@
 !
 !-- Flag parameters for RRTMGS (should not be changed)
-    INTEGER(iwp), PARAMETER :: rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
+    INTEGER(iwp), PARAMETER :: rrtm_idrv     = 1, & !< flag for longwave upward flux calculation option (0,1)
+                               rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
                                rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
                                rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
@@ -289,6 +292,5 @@
     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 = 1        !< longwave upward flux calculation option (0,1)
+                    rrtm_iaer = 0        !< aerosol option flag (0: no aerosol layers, for lw only: 6 (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented)
 
     INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
@@ -385,4 +387,28 @@
     END INTERFACE radiation_constant
  
+    INTERFACE radiation_control
+       MODULE PROCEDURE radiation_control
+    END INTERFACE radiation_control
+
+    INTERFACE radiation_3d_data_averaging
+       MODULE PROCEDURE radiation_3d_data_averaging
+    END INTERFACE radiation_3d_data_averaging
+
+    INTERFACE radiation_data_output_2d
+       MODULE PROCEDURE radiation_data_output_2d
+    END INTERFACE radiation_data_output_2d
+
+    INTERFACE radiation_data_output_3d
+       MODULE PROCEDURE radiation_data_output_3d
+    END INTERFACE radiation_data_output_3d
+
+    INTERFACE radiation_data_output_mask
+       MODULE PROCEDURE radiation_data_output_mask
+    END INTERFACE radiation_data_output_mask
+
+    INTERFACE radiation_define_netcdf_grid
+       MODULE PROCEDURE radiation_define_netcdf_grid
+    END INTERFACE radiation_define_netcdf_grid
+
     INTERFACE radiation_header
        MODULE PROCEDURE radiation_header
@@ -406,4 +432,12 @@
     END INTERFACE radiation_tendency
 
+    INTERFACE radiation_read_restart_data
+       MODULE PROCEDURE radiation_read_restart_data
+    END INTERFACE radiation_read_restart_data
+
+    INTERFACE radiation_last_actions
+       MODULE PROCEDURE radiation_last_actions
+    END INTERFACE radiation_last_actions
+
     SAVE
 
@@ -411,12 +445,15 @@
 
 !
-!-- Public functions
+!-- Public functions / NEEDS SORTING
     PUBLIC radiation_check_data_output, radiation_check_data_output_pr,        &
-           radiation_check_parameters, radiation_clearsky, radiation_constant, &
-           radiation_header, radiation_init, radiation_parin, radiation_rrtmg, &
-           radiation_tendency
+           radiation_check_parameters, radiation_control,                      &
+           radiation_header, radiation_init, radiation_parin,                  &
+           radiation_3d_data_averaging, radiation_tendency,                    &
+           radiation_data_output_2d, radiation_data_output_3d,                 &
+           radiation_define_netcdf_grid, radiation_last_actions,               &
+           radiation_read_restart_data, radiation_data_output_mask
     
 !
-!-- Public variables and constants
+!-- Public variables and constants / NEEDS SORTING
     PUBLIC dots_rad, dt_radiation, force_radiation_call,                       &
            rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
@@ -424,13 +461,43 @@
            rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
            rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
-           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb,                 &
+           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av,                           &
            skip_time_do_radiation, time_radiation, unscheduled_radiation_calls
 
 
 #if defined ( __rrtmg )
-    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv
+    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
 #endif
 
  CONTAINS
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine controls the calls of the radiation schemes
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_control
+ 
+ 
+       IMPLICIT NONE
+
+
+       SELECT CASE ( TRIM( radiation_scheme ) )
+
+          CASE ( 'constant' )
+             CALL radiation_constant
+          
+          CASE ( 'clear-sky' )
+             CALL radiation_clearsky
+       
+          CASE ( 'rrtmg' )
+             CALL radiation_rrtmg
+
+          CASE DEFAULT
+
+       END SELECT
+
+
+    END SUBROUTINE radiation_control
 
 !------------------------------------------------------------------------------!
@@ -456,6 +523,5 @@
        SELECT CASE ( TRIM( var ) )
 
-         CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_cs_hr', 'rad_lw_hr',       &
-                 'rad_sw_in', 'rad_sw_out', 'rad_sw_cs_hr', 'rad_sw_hr' )
+         CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr' )
              IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
                 message_string = '"output of "' // TRIM( var ) // '" requi' // &
@@ -779,14 +845,6 @@
        ENDIF
 
-
-       IF ( radiation_scheme == 'constant' )  THEN
-
-          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
-             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
-          ENDIF
-
-       ENDIF
-
-       IF ( radiation_scheme == 'clear-sky' )  THEN
+       IF ( radiation_scheme == 'clear-sky'  .OR.                              &
+            radiation_scheme == 'constant')  THEN
 
           ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
@@ -1075,4 +1133,8 @@
                             + rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
 
+
+             rad_lw_out_change_0(j,i) = 3.0_wp * sigma_sb * emissivity         &
+                                        * (pt(k,j,i) * exn) ** 3
+
           ENDDO
        ENDDO
@@ -1093,4 +1155,5 @@
        INTEGER(iwp) :: i, j, k   !< loop indices
        REAL(wp)     :: exn,   &  !< Exner functions at surface
+                       exn1,  &  !< Exner functions at first grid level
                        pt1       !< potential temperature at first grid level
 
@@ -1099,6 +1162,5 @@
        exn = (surface_pressure / 1000.0_wp )**0.286_wp
 !
-!--    Prescribe net radiation and estimate a longwave outgoing radiative
-!--    flux (needed in land surface model)
+!--    Prescribe net radiation and estimate the remaining radiative fluxes
        DO i = nxlg, nxrg
           DO j = nysg, nyng
@@ -1106,5 +1168,19 @@
 
              rad_net(j,i)      = net_radiation
+
+             exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
+
+             IF ( cloud_physics )  THEN
+                pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
+                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
+             ELSE
+                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
+             ENDIF
+
              rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
+
+             rad_sw_in(0,j,i) = ( rad_net(j,i) - rad_lw_in(0,j,i)              &
+                                  + rad_lw_out(0,j,i) )                        &
+                                  / ( 1.0_wp - alpha(j,i) )
 
           ENDDO
@@ -2147,24 +2223,25 @@
 !> Cache-optimized version.
 !------------------------------------------------------------------------------!
-    SUBROUTINE radiation_tendency_ij ( i, j, tend )
-
-       USE cloud_parameters,                                                   &
-           ONLY:  pt_d_t
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) :: i, j, k !< loop indices
-
-       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
-
-#if defined ( __rrtmg )
+ SUBROUTINE radiation_tendency_ij ( i, j, tend )
+
+    USE cloud_parameters,                                                      &
+        ONLY:  pt_d_t
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) :: i, j, k !< loop indices
+
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
+
+    IF ( radiation_scheme == 'rrtmg' )  THEN
+#if defined  ( __rrtmg )
 !
 !--    Calculate tendency based on heating rate
        DO k = nzb+1, nzt+1
           tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i))    &
-                                      * pt_d_t(k) * d_seconds_hour
+                                         * pt_d_t(k) * d_seconds_hour
        ENDDO
-
 #endif
+    ENDIF
 
     END SUBROUTINE radiation_tendency_ij
@@ -2177,19 +2254,20 @@
 !> Vector-optimized version
 !------------------------------------------------------------------------------!
-    SUBROUTINE radiation_tendency ( tend )
-
-       USE cloud_parameters,                                                   &
-           ONLY:  pt_d_t
-
-       USE indices,                                                            &
-           ONLY:  nxl, nxr, nyn, nys
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) :: i, j, k !< loop indices
-
-       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
-
-#if defined ( __rrtmg )
+ SUBROUTINE radiation_tendency ( tend )
+
+    USE cloud_parameters,                                                      &
+        ONLY:  pt_d_t
+
+    USE indices,                                                               &
+        ONLY:  nxl, nxr, nyn, nys
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) :: i, j, k !< loop indices
+
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
+
+    IF ( radiation_scheme == 'rrtmg' )  THEN
+#if defined  ( __rrtmg )
 !
 !--    Calculate tendency based on heating rate
@@ -2198,12 +2276,1280 @@
              DO k = nzb+1, nzt+1
                 tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i)                 &
-                                            +  rad_sw_hr(k,j,i) ) * pt_d_t(k)  &
-                                            * d_seconds_hour
-             ENDDO
-         ENDDO
+                                          +  rad_sw_hr(k,j,i) ) * pt_d_t(k)    &
+                                          * d_seconds_hour
+             ENDDO
+          ENDDO
        ENDDO
 #endif
-
-    END SUBROUTINE radiation_tendency
+    ENDIF
+
+
+ END SUBROUTINE radiation_tendency
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine for averaging 3D data
+!------------------------------------------------------------------------------!
+SUBROUTINE radiation_3d_data_averaging( mode, variable )
+ 
+
+    USE control_parameters
+
+    USE indices
+
+    USE kinds
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  mode    !< 
+    CHARACTER (LEN=*) :: variable !< 
+
+    INTEGER(iwp) ::  i !< 
+    INTEGER(iwp) ::  j !< 
+    INTEGER(iwp) ::  k !< 
+
+    IF ( mode == 'allocate' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+             CASE ( 'rad_net*' )
+                IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
+                   ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_net_av = 0.0_wp
+
+             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_lw_cs_hr' )
+                IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
+                   ALLOCATE( rad_lw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_lw_cs_hr_av = 0.0_wp
+
+             CASE ( 'rad_lw_hr' )
+                IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
+                   ALLOCATE( rad_lw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_lw_hr_av = 0.0_wp
+
+             CASE ( 'rad_sw_in' )
+                IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
+                   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 ( 'rad_sw_cs_hr' )
+                IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
+                   ALLOCATE( rad_sw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_sw_cs_hr_av = 0.0_wp
+
+             CASE ( 'rad_sw_hr' )
+                IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
+                   ALLOCATE( rad_sw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_sw_hr_av = 0.0_wp
+
+          CASE DEFAULT
+             CONTINUE
+
+       END SELECT
+
+    ELSEIF ( mode == 'sum' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+          CASE ( 'rad_net*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   rad_net_av(j,i) = rad_net_av(j,i) + rad_net(j,i)
+                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) + 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_lw_cs_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) + rad_lw_cs_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_lw_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) + rad_lw_hr(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
+
+          CASE ( 'rad_sw_cs_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) + rad_sw_cs_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) + rad_sw_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE DEFAULT
+             CONTINUE
+
+       END SELECT
+
+    ELSEIF ( mode == 'average' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+         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_lw_cs_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_lw_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_hr_av(k,j,i) = rad_lw_hr_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
+
+          CASE ( 'rad_sw_cs_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+       END SELECT
+
+    ENDIF
+
+END SUBROUTINE radiation_3d_data_averaging
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining appropriate grid for netcdf variables.
+!> It is called out from subroutine netcdf.
+!------------------------------------------------------------------------------!
+SUBROUTINE radiation_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
+    
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*), INTENT(IN)  ::  var         !< 
+    LOGICAL, INTENT(OUT)           ::  found       !< 
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !< 
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !< 
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !< 
+
+    found  = .TRUE.
+
+
+!
+!-- Check for the grid
+    SELECT CASE ( TRIM( var ) )
+
+       CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr',        &
+              'rad_lw_cs_hr_xy', 'rad_lw_hr_xy', 'rad_sw_cs_hr_xy',            &
+              'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',               &
+              'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz',            &
+              'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz' )
+          grid_x = 'x'
+          grid_y = 'y'
+          grid_z = 'zu'
+
+       CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',            &
+              'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy','rad_sw_out_xy', &
+              'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz','rad_sw_out_xz', &
+              'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz','rad_sw_out_yz' )
+          grid_x = 'x'
+          grid_y = 'y'
+          grid_z = 'zw'
+
+
+       CASE DEFAULT
+          found  = .FALSE.
+          grid_x = 'none'
+          grid_y = 'none'
+          grid_z = 'none'
+
+        END SELECT
+
+    END SUBROUTINE radiation_define_netcdf_grid
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining 3D output variables
+!------------------------------------------------------------------------------!
+ SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode,         &
+                                      local_pf, two_d )
+ 
+    USE indices
+
+    USE kinds
+
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  grid     !< 
+    CHARACTER (LEN=*) ::  mode     !< 
+    CHARACTER (LEN=*) ::  variable !< 
+
+    INTEGER(iwp) ::  av !< 
+    INTEGER(iwp) ::  i  !< 
+    INTEGER(iwp) ::  j  !< 
+    INTEGER(iwp) ::  k  !< 
+
+    LOGICAL      ::  found !< 
+    LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
+
+    REAL(wp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !< 
+
+    found = .TRUE.
+
+    SELECT CASE ( TRIM( variable ) )
+
+       CASE ( 'rad_net*_xy' )        ! 2d-array
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = rad_net(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   local_pf(i,j,nzb+1) = rad_net_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+          two_d = .TRUE.
+          grid = 'zu1'
+
+ 
+       CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_in(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_in_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zu'
+
+       CASE ( 'rad_lw_out_xy', 'rad_lw_out_xz', 'rad_lw_out_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_out(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_out_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF   
+          IF ( mode == 'xy' )  grid = 'zu'
+
+       CASE ( 'rad_lw_cs_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_cs_hr_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zw'
+
+       CASE ( 'rad_lw_hr_xy', 'rad_lw_hr_xz', 'rad_lw_hr_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zw'
+
+       CASE ( 'rad_sw_in_xy', 'rad_sw_in_xz', 'rad_sw_in_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_in(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_in_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zu'
+
+       CASE ( 'rad_sw_out_xy', 'rad_sw_out_xz', 'rad_sw_out_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_out(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_out_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zu'
+
+       CASE ( 'rad_sw_cs_hr_xy', 'rad_sw_cs_hr_xz', 'rad_sw_cs_hr_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zw'
+
+       CASE ( 'rad_sw_hr_xy', 'rad_sw_hr_xz', 'rad_sw_hr_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zw'
+
+       CASE DEFAULT
+          found = .FALSE.
+          grid  = 'none'
+
+    END SELECT
+ 
+ END SUBROUTINE radiation_data_output_2d
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining 3D output variables
+!------------------------------------------------------------------------------!
+ SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf )
+ 
+
+    USE indices
+
+    USE kinds
+
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  variable !< 
+
+    INTEGER(iwp) ::  av    !< 
+    INTEGER(iwp) ::  i     !< 
+    INTEGER(iwp) ::  j     !< 
+    INTEGER(iwp) ::  k     !< 
+
+    LOGICAL      ::  found !< 
+
+    REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !< 
+
+
+    found = .TRUE.
+
+
+    SELECT CASE ( TRIM( variable ) )
+
+      CASE ( 'rad_sw_in' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_in(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_in_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_sw_out' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_out(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_out_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_sw_cs_hr' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_sw_hr' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_hr(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_lw_in' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_in(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_in_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_lw_out' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_out(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_out_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_lw_cs_hr' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_lw_hr' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_hr(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+       CASE DEFAULT
+          found = .FALSE.
+
+    END SELECT
+
+
+ END SUBROUTINE radiation_data_output_3d
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining masked data output
+!------------------------------------------------------------------------------!
+ SUBROUTINE radiation_data_output_mask( av, variable, found, local_pf )
+ 
+    USE control_parameters
+        
+    USE indices
+    
+    USE kinds
+    
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  variable   !< 
+
+    INTEGER(iwp) ::  av   !< 
+    INTEGER(iwp) ::  i    !< 
+    INTEGER(iwp) ::  j    !< 
+    INTEGER(iwp) ::  k    !< 
+
+    LOGICAL ::  found     !< 
+
+    REAL(wp),                                                                  &
+       DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
+          local_pf   !< 
+
+
+    found = .TRUE.
+
+    SELECT CASE ( TRIM( variable ) )
+
+
+       CASE ( 'rad_lw_in' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_in(mask_k(mid,k),              &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_in_av(mask_k(mid,k),           &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_lw_out' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_out(mask_k(mid,k),             &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_out_av(mask_k(mid,k),          &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_lw_cs_hr' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_cs_hr(mask_k(mid,k),           &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_cs_hr_av(mask_k(mid,k),        &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_lw_hr' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_hr(mask_k(mid,k),              &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_hr_av(mask_k(mid,k),           &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_sw_in' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_in(mask_k(mid,k),              &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_in_av(mask_k(mid,k),           &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_sw_out' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_out(mask_k(mid,k),             &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_out_av(mask_k(mid,k),          &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_sw_cs_hr' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_cs_hr(mask_k(mid,k),           &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_cs_hr_av(mask_k(mid,k),        &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_sw_hr' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_hr(mask_k(mid,k),              &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_hr_av(mask_k(mid,k),           &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE DEFAULT
+          found = .FALSE.
+
+    END SELECT
+
+
+ END SUBROUTINE radiation_data_output_mask
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defines masked output variables
+!------------------------------------------------------------------------------!
+ SUBROUTINE radiation_last_actions
+ 
+
+    USE control_parameters
+        
+    USE kinds
+
+    IMPLICIT NONE
+
+    IF ( write_binary(1:4) == 'true' )  THEN
+       IF ( ALLOCATED( rad_net ) )  THEN
+          WRITE ( 14 )  'rad_net             ';  WRITE ( 14 )  rad_net  
+       ENDIF
+       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_lw_out_change_0 ) )  THEN
+          WRITE ( 14 )  'rad_lw_out_change_0 '
+          WRITE ( 14 )  rad_lw_out_change_0
+       ENDIF
+       IF ( ALLOCATED( rad_lw_cs_hr ) )  THEN
+          WRITE ( 14 )  'rad_lw_cs_hr        ';  WRITE ( 14 )  rad_lw_cs_hr
+       ENDIF
+       IF ( ALLOCATED( rad_lw_cs_hr_av ) )  THEN
+          WRITE ( 14 )  'rad_lw_cs_hr_av     ';  WRITE ( 14 )  rad_lw_cs_hr_av
+       ENDIF
+       IF ( ALLOCATED( rad_lw_hr ) )  THEN
+          WRITE ( 14 )  'rad_lw_hr           ';  WRITE ( 14 )  rad_lw_hr
+       ENDIF
+       IF ( ALLOCATED( rad_lw_hr_av ) )  THEN
+          WRITE ( 14 )  'rad_lw_hr_av        ';  WRITE ( 14 )  rad_lw_hr_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  
+       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 
+       IF ( ALLOCATED( rad_sw_cs_hr ) )  THEN
+          WRITE ( 14 )  'rad_sw_cs_hr        ';  WRITE ( 14 )  rad_sw_cs_hr
+       ENDIF
+       IF ( ALLOCATED( rad_sw_cs_hr_av ) )  THEN
+          WRITE ( 14 )  'rad_sw_cs_hr_av     ';  WRITE ( 14 )  rad_sw_cs_hr_av
+       ENDIF
+       IF ( ALLOCATED( rad_sw_hr ) )  THEN
+          WRITE ( 14 )  'rad_sw_hr           ';  WRITE ( 14 )  rad_sw_hr
+       ENDIF
+       IF ( ALLOCATED( rad_sw_hr_av ) )  THEN
+          WRITE ( 14 )  'rad_sw_hr_av        ';  WRITE ( 14 )  rad_sw_hr_av
+       ENDIF
+
+       WRITE ( 14 )  '*** end rad ***     '
+
+    ENDIF
+
+ END SUBROUTINE radiation_last_actions
+
+
+SUBROUTINE radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file,   &
+                                     nynfa, nyn_on_file, nysfa, nys_on_file,   &
+                                     offset_xa, offset_ya, overlap_count,      &
+                                     tmp_2d, tmp_3d )
+ 
+
+    USE control_parameters
+        
+    USE indices
+    
+    USE kinds
+    
+    USE pegrid
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=20) :: field_char   !< 
+
+    INTEGER(iwp) ::  i               !< 
+    INTEGER(iwp) ::  k               !< 
+    INTEGER(iwp) ::  nxlc            !< 
+    INTEGER(iwp) ::  nxlf            !< 
+    INTEGER(iwp) ::  nxl_on_file     !< 
+    INTEGER(iwp) ::  nxrc            !< 
+    INTEGER(iwp) ::  nxrf            !< 
+    INTEGER(iwp) ::  nxr_on_file     !< 
+    INTEGER(iwp) ::  nync            !< 
+    INTEGER(iwp) ::  nynf            !< 
+    INTEGER(iwp) ::  nyn_on_file     !< 
+    INTEGER(iwp) ::  nysc            !< 
+    INTEGER(iwp) ::  nysf            !< 
+    INTEGER(iwp) ::  nys_on_file     !< 
+    INTEGER(iwp) ::  overlap_count   !< 
+
+    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxlfa       !< 
+    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxrfa       !< 
+    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nynfa       !< 
+    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nysfa       !< 
+    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_xa   !< 
+    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_ya   !< 
+
+    REAL(wp),                                                                  &
+       DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
+          tmp_2d   !< 
+
+    REAL(wp),                                                                  &
+       DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
+          tmp_3d   !< 
+
+
+
+   IF ( initializing_actions == 'read_restart_data' )  THEN
+      READ ( 13 )  field_char
+
+      DO  WHILE ( TRIM( field_char ) /= '*** end rad ***' )
+
+         DO  k = 1, overlap_count
+
+            nxlf = nxlfa(i,k)
+            nxlc = nxlfa(i,k) + offset_xa(i,k)
+            nxrf = nxrfa(i,k)
+            nxrc = nxrfa(i,k) + offset_xa(i,k)
+            nysf = nysfa(i,k)
+            nysc = nysfa(i,k) + offset_ya(i,k)
+            nynf = nynfa(i,k)
+            nync = nynfa(i,k) + offset_ya(i,k)
+
+
+            SELECT CASE ( TRIM( field_char ) )
+
+                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
+                      ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                CASE ( 'rad_lw_in' )
+                   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_lw_out_change_0' )
+                   IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) )  THEN
+                      ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)&
+                              = tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_cs_hr' )
+                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr ) )  THEN
+                      ALLOCATE( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_cs_hr_av' )
+                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
+                      ALLOCATE( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_hr' )
+                   IF ( .NOT. ALLOCATED( rad_lw_hr ) )  THEN
+                      ALLOCATE( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_hr_av' )
+                   IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
+                      ALLOCATE( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_hr_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(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 ( 'rad_sw_cs_hr' )
+                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr ) )  THEN
+                      ALLOCATE( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_cs_hr_av' )
+                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
+                      ALLOCATE( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_hr' )
+                   IF ( .NOT. ALLOCATED( rad_sw_hr ) )  THEN
+                      ALLOCATE( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_hr_av' )
+                   IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
+                      ALLOCATE( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+               CASE DEFAULT
+                  WRITE( message_string, * ) 'unknown variable named "',       &
+                                        TRIM( field_char ), '" found in',      &
+                                        '&data from prior run on PE ', myid
+                  CALL message( 'radiation_read_restart_data', 'PA0441', 1, 2, 0, 6, &
+                                 0 )
+
+            END SELECT
+
+         ENDDO
+
+         READ ( 13 )  field_char
+
+      ENDDO
+   ENDIF
+
+ END SUBROUTINE radiation_read_restart_data
+
 
  END MODULE radiation_model_mod
Index: palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/read_3d_binary.f90	(revision 1975)
+++ palm/trunk/SOURCE/read_3d_binary.f90	(revision 1976)
@@ -19,5 +19,7 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix: read of land surface data only when module is switched on.
+! Radiation parts are now done in the respective module.
+! Binary version increased to 4.5.
 ! 
 ! Former revisions:
@@ -138,9 +140,5 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY: rad_net, rad_net_av, rad_lw_in, rad_lw_in_av, rad_lw_out,        &
-              rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av,          &
-              rad_lw_out_av, rad_lw_out_change_0, rad_sw_in, rad_sw_in_av,     &
-              rad_sw_out, rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av,        &
-              rad_sw_hr, rad_sw_hr_av
+        ONLY: radiation, radiation_read_restart_data
 
     USE random_function_mod,                                                   &
@@ -328,5 +326,5 @@
 !--    First compare the version numbers
        READ ( 13 )  version_on_file
-       binary_version = '4.4'
+       binary_version = '4.5'
        IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
           WRITE( message_string, * ) 'version mismatch concerning data ',      &
@@ -698,155 +696,4 @@
                    qv_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
                              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
-                      ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF  
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
-                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-                CASE ( 'rad_lw_in' )
-                   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_lw_out_change_0' )
-                   IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) )  THEN
-                      ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)&
-                              = tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'rad_lw_cs_hr' )
-                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr ) )  THEN
-                      ALLOCATE( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
-                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'rad_lw_cs_hr_av' )
-                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
-                      ALLOCATE( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_lw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
-                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'rad_lw_hr' )
-                   IF ( .NOT. ALLOCATED( rad_lw_hr ) )  THEN
-                      ALLOCATE( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
-                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'rad_lw_hr_av' )
-                   IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
-                      ALLOCATE( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_lw_hr_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(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 ( 'rad_sw_cs_hr' )
-                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr ) )  THEN
-                      ALLOCATE( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
-                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'rad_sw_cs_hr_av' )
-                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
-                      ALLOCATE( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_sw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
-                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'rad_sw_hr' )
-                   IF ( .NOT. ALLOCATED( rad_sw_hr ) )  THEN
-                      ALLOCATE( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
-                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'rad_sw_hr_av' )
-                   IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
-                      ALLOCATE( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_lw_hr_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
@@ -1274,8 +1121,20 @@
 !
 !--    Read land surface restart data
-       CALL lsm_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file,  &
-                                    nynfa, nyn_on_file, nysfa, nys_on_file,    &
-                                    offset_xa, offset_ya, overlap_count(i),    &
-                                    tmp_2d )
+       IF ( land_surface )  THEN
+          CALL lsm_read_restart_data( i, nxlfa, nxl_on_file, nxrfa,            &
+                                      nxr_on_file, nynfa, nyn_on_file, nysfa,  &
+                                      nys_on_file, offset_xa, offset_ya,       &
+                                      overlap_count(i), tmp_2d )
+       ENDIF
+
+!
+!--    Read radiation restart data
+       IF ( radiation )  THEN
+          CALL radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa,      &
+                                            nxr_on_file, nynfa, nyn_on_file,   &
+                                            nysfa, nys_on_file, offset_xa,     &
+                                            offset_ya, overlap_count(i),       &
+                                            tmp_2d, tmp_3d )
+       ENDIF
 
 !
Index: palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 1975)
+++ palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 1976)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Radiation actions are now done directly in the respective module
 ! 
 ! Former revisions:
@@ -138,9 +138,5 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av, rad_sw_out,       &
-               rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,        &
-               rad_sw_hr_av, rad_lw_in, rad_lw_in_av, rad_lw_out,              &
-               rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr,        &
-               rad_lw_hr_av
+        ONLY:  radiation, radiation_3d_data_averaging
 
 
@@ -291,58 +287,4 @@
                 qv_av = 0.0_wp
 
-             CASE ( 'rad_net*' )
-                IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
-                   ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                rad_net_av = 0.0_wp
-
-             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_lw_cs_hr' )
-                IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
-                   ALLOCATE( rad_lw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                rad_lw_cs_hr_av = 0.0_wp
-
-             CASE ( 'rad_lw_hr' )
-                IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
-                   ALLOCATE( rad_lw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                rad_lw_hr_av = 0.0_wp
-
-             CASE ( 'rad_sw_in' )
-                IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
-                   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 ( 'rad_sw_cs_hr' )
-                IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
-                   ALLOCATE( rad_sw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                rad_sw_cs_hr_av = 0.0_wp
-
-             CASE ( 'rad_sw_hr' )
-                IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
-                   ALLOCATE( rad_sw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                rad_sw_hr_av = 0.0_wp
-
              CASE ( 'rho' )
                 IF ( .NOT. ALLOCATED( rho_av ) )  THEN
@@ -429,4 +371,10 @@
                 IF ( land_surface )  THEN
                    CALL lsm_3d_data_averaging( 'allocate', doav(ii) )
+                ENDIF
+
+!
+!--             Radiation quantity
+                IF ( radiation )  THEN
+                   CALL radiation_3d_data_averaging( 'allocate', doav(ii) )
                 ENDIF
 
@@ -655,83 +603,4 @@
              ENDDO
 
-          CASE ( 'rad_net*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   rad_net_av(j,i) = rad_net_av(j,i) + rad_net(j,i)
-                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) + 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_lw_cs_hr' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) + rad_lw_cs_hr(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
-
-          CASE ( 'rad_lw_hr' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) + rad_lw_hr(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
-
-          CASE ( 'rad_sw_cs_hr' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) + rad_sw_cs_hr(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
-
-          CASE ( 'rad_sw_hr' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
-                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) + rad_sw_hr(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
-
           CASE ( 'rho' )
              DO  i = nxlg, nxrg
@@ -854,4 +723,10 @@
 
 !
+!--          Radiation quantity
+             IF ( radiation )  THEN
+                CALL radiation_3d_data_averaging( 'sum', doav(ii) )
+             ENDIF
+
+!
 !--          User-defined quantity
              CALL user_3d_data_averaging( 'sum', doav(ii) )
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 1975)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 1976)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Simplified calls to radiation model
 ! 
 ! Former revisions:
@@ -321,7 +321,6 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY: dt_radiation, force_radiation_call, radiation,                   &
-              radiation_clearsky, radiation_constant, radiation_rrtmg,         &
-              radiation_scheme, skip_time_do_radiation, time_radiation
+        ONLY: dt_radiation, force_radiation_call, radiation, radiation_control,&
+              skip_time_do_radiation, time_radiation
 
     USE spectra_mod,                                                           &
@@ -900,11 +899,5 @@
                 ENDIF
 
-                IF ( radiation_scheme == 'clear-sky' )  THEN
-                   CALL radiation_clearsky
-                ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
-                   CALL radiation_rrtmg
-                ELSE
-                   CALL radiation_constant
-                ENDIF
+                CALL radiation_control
 
                 CALL cpu_log( log_point(50), 'radiation', 'stop' )
Index: palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/write_3d_binary.f90	(revision 1975)
+++ palm/trunk/SOURCE/write_3d_binary.f90	(revision 1976)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Radiation actions are now done directly in the respective module.
+! Binary version increased to 4.5.
 ! 
 ! Former revisions:
@@ -122,11 +123,4 @@
     USE pegrid
     
-    USE radiation_model_mod,                                                   &
-        ONLY: radiation, rad_net, rad_net_av, rad_lw_in, rad_lw_in_av,         &
-              rad_lw_out, rad_lw_out_av, rad_lw_out_change_0, rad_lw_cs_hr,   &
-              rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,             &
-              rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,           &
-              rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av
-
     USE random_function_mod,                                                   &
         ONLY:  random_iv, random_iy
@@ -149,5 +143,5 @@
 !
 !-- Write arrays.
-    binary_version = '4.4'
+    binary_version = '4.5'
 
     WRITE ( 14 )  binary_version
@@ -254,64 +248,4 @@
        WRITE ( 14 )  'sswst               ';  WRITE ( 14 ) sswst
     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_lw_out_change_0 ) )  THEN
-          WRITE ( 14 )  'rad_lw_out_change_0 '
-          WRITE ( 14 )  rad_lw_out_change_0
-       ENDIF
-       IF ( ALLOCATED( rad_lw_cs_hr ) )  THEN
-          WRITE ( 14 )  'rad_lw_cs_hr        ';  WRITE ( 14 )  rad_lw_cs_hr
-       ENDIF
-       IF ( ALLOCATED( rad_lw_cs_hr_av ) )  THEN
-          WRITE ( 14 )  'rad_lw_cs_hr_av     ';  WRITE ( 14 )  rad_lw_cs_hr_av
-       ENDIF
-       IF ( ALLOCATED( rad_lw_hr ) )  THEN
-          WRITE ( 14 )  'rad_lw_hr           ';  WRITE ( 14 )  rad_lw_hr
-       ENDIF
-       IF ( ALLOCATED( rad_lw_hr_av ) )  THEN
-          WRITE ( 14 )  'rad_lw_hr_av        ';  WRITE ( 14 )  rad_lw_hr_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  
-       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 
-       IF ( ALLOCATED( rad_sw_cs_hr ) )  THEN
-          WRITE ( 14 )  'rad_sw_cs_hr        ';  WRITE ( 14 )  rad_sw_cs_hr
-       ENDIF
-       IF ( ALLOCATED( rad_sw_cs_hr_av ) )  THEN
-          WRITE ( 14 )  'rad_sw_cs_hr_av     ';  WRITE ( 14 )  rad_sw_cs_hr_av
-       ENDIF
-       IF ( ALLOCATED( rad_sw_hr ) )  THEN
-          WRITE ( 14 )  'rad_sw_hr           ';  WRITE ( 14 )  rad_sw_hr
-       ENDIF
-       IF ( ALLOCATED( rad_sw_hr_av ) )  THEN
-          WRITE ( 14 )  'rad_sw_hr_av        ';  WRITE ( 14 )  rad_sw_hr_av
-       ENDIF
-    ENDIF
     IF ( ocean )  THEN
        IF ( ALLOCATED( rho_av ) )  THEN
