Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 1971)
+++ /palm/trunk/SOURCE/Makefile	(revision 1972)
@@ -20,5 +20,5 @@
 # Current revisions:
 # ------------------
-# 
+# Removed some dependencies due to further modularization of land surface model
 # 
 # Former revisions:
@@ -509,5 +509,5 @@
 nudging_mod.o: modules.o cpulog_mod.o mod_kinds.o
 package_parin.o: modules.o mod_kinds.o mod_particle_attributes.o
-palm.o: modules.o cpulog_mod.o ls_forcing_mod.o mod_kinds.o nudging_mod.o\
+palm.o: modules.o cpulog_mod.o land_surface_model_mod.o ls_forcing_mod.o mod_kinds.o nudging_mod.o\
 	pmc_interface_mod.o surface_layer_fluxes_mod.o
 parin.o: modules.o cpulog_mod.o land_surface_model_mod.o mod_kinds.o netcdf_interface_mod.o \
@@ -541,6 +541,6 @@
 random_gauss.o: mod_kinds.o random_function_mod.o random_generator_parallel_mod.o
 random_generator_parallel_mod.o: mod_kinds.o
-read_3d_binary.o: modules.o cpulog_mod.o mod_kinds.o land_surface_model_mod.o \
-   radiation_model_mod.o random_function_mod.o random_generator_parallel_mod.o \
+read_3d_binary.o: modules.o cpulog_mod.o mod_kinds.o \
+   land_surface_model_mod.o radiation_model_mod.o random_function_mod.o random_generator_parallel_mod.o \
    spectra_mod.o
 read_var_list.o: modules.o mod_kinds.o netcdf_interface_mod.o plant_canopy_model_mod.o \
@@ -604,5 +604,5 @@
 wall_fluxes.o: modules.o mod_kinds.o
 wind_turbine_model_mod.o: modules.o cpulog_mod.o mod_kinds.o
-write_3d_binary.o: modules.o cpulog_mod.o mod_kinds.o land_surface_model_mod.o \
+write_3d_binary.o: modules.o cpulog_mod.o mod_kinds.o \
 	radiation_model_mod.o random_function_mod.o random_generator_parallel_mod.o \
 	spectra_mod.o
Index: /palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- /palm/trunk/SOURCE/average_3d_data.f90	(revision 1971)
+++ /palm/trunk/SOURCE/average_3d_data.f90	(revision 1972)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Output of land surface quantities is now done directly in the respective module
 ! 
 ! Former revisions:
@@ -95,8 +95,5 @@
 
     USE land_surface_model_mod,                                                &
-        ONLY:  c_liq_av, c_soil_av, c_veg_av, ghf_eb_av, lai_av, m_liq_eb_av,  &
-               m_soil_av, nzb_soil, nzt_soil, qsws_eb_av, qsws_liq_eb_av,      &
-               qsws_soil_eb_av, qsws_veg_eb_av, r_a, r_a_av, r_s, r_s_av,      &
-               shf_eb_av, t_soil_av
+        ONLY:  land_surface, lsm_3d_data_averaging
 
     USE radiation_model_mod,                                                   &
@@ -130,25 +127,4 @@
        SELECT CASE ( TRIM( doav(ii) ) )
 
-         CASE ( 'c_liq*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   c_liq_av(j,i) = c_liq_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
-         CASE ( 'c_soil*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   c_soil_av(j,i) = c_soil_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
-         CASE ( 'c_veg*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   c_veg_av(j,i) = c_veg_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
           CASE ( 'e' )
              DO  i = nxlg, nxrg
@@ -160,11 +136,4 @@
              ENDDO
 
-         CASE ( 'ghf_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   ghf_eb_av(j,i) = ghf_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
           CASE ( 'qsws*' )
              DO  i = nxlg, nxrg
@@ -174,11 +143,4 @@
              ENDDO
 
-         CASE ( 'lai*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   lai_av(j,i) = lai_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
           CASE ( 'lpt' )
              DO  i = nxlg, nxrg
@@ -194,20 +156,4 @@
                 DO  j = nysg, nyng
                    lwp_av(j,i) = lwp_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
-         CASE ( 'm_liq_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   m_liq_eb_av(j,i) = m_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
-          CASE ( 'm_soil' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb_soil, nzt_soil
-                      m_soil_av(k,j,i) = m_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
-                   ENDDO
                 ENDDO
              ENDDO
@@ -337,32 +283,4 @@
              ENDDO
 
-         CASE ( 'qsws_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   qsws_eb_av(j,i) = qsws_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
-         CASE ( 'qsws_liq_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
-         CASE ( 'qsws_soil_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
-         CASE ( 'qsws_veg_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
           CASE ( 'qv' )
              DO  i = nxlg, nxrg
@@ -453,18 +371,4 @@
              ENDDO
 
-         CASE ( 'r_a*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   r_a_av(j,i) = r_a_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
-         CASE ( 'r_s*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   r_s_av(j,i) = r_s_av(j,i) / REAL( average_count_3d, KIND=wp )
-                ENDDO
-             ENDDO
-
           CASE ( 'rho' )
              DO  i = nxlg, nxrg
@@ -515,13 +419,4 @@
              ENDDO
 
-          CASE ( 't_soil' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb_soil, nzt_soil
-                      t_soil_av(k,j,i) = t_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
-                   ENDDO
-                ENDDO
-             ENDDO
-
           CASE ( 'u' )
              DO  i = nxlg, nxrg
@@ -582,4 +477,10 @@
 
           CASE DEFAULT
+!
+!--          Land surface quantity
+             IF ( land_surface )  THEN
+                CALL lsm_3d_data_averaging( 'average', doav(ii) )
+             ENDIF
+
 !
 !--          User-defined quantity
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 1971)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 1972)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Removed check of lai* (done in land surface module)
 ! 
 ! Former revisions:
@@ -2960,5 +2960,5 @@
              CONTINUE
 
-          CASE ( 'lai*', 'lwp*', 'ol*', 'pra*', 'prr*', 'qsws*', 'shf*', 't*', &
+          CASE ( 'lwp*', 'ol*', 'pra*', 'prr*', 'qsws*', 'shf*', 't*', &
                  'u*', 'z0*', 'z0h*', 'z0q*' )
              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
Index: /palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_2d.f90	(revision 1971)
+++ /palm/trunk/SOURCE/data_output_2d.f90	(revision 1972)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Output of land surface quantities is now done directly in the respective module
 ! 
 ! Former revisions:
@@ -174,10 +174,5 @@
     
     USE land_surface_model_mod,                                                &
-        ONLY:  c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb,            &
-               ghf_eb_av, lai, lai_av, m_liq_eb, m_liq_eb_av, m_soil,          &
-               m_soil_av, nzb_soil, nzt_soil, qsws_eb, qsws_eb_av,             &
-               qsws_liq_eb, qsws_liq_eb_av, qsws_soil_eb, qsws_soil_eb_av,     &
-               qsws_veg_eb, qsws_veg_eb_av, r_a, r_a_av, r_s, r_s_av, shf_eb,  &
-               shf_eb_av, t_soil, t_soil_av, zs
+        ONLY:  land_surface, lsm_data_output_2d, zs
     
 #if defined( __netcdf )
@@ -290,5 +285,4 @@
              CALL check_open( 101+av*10 )
           ENDIF
-
           IF ( data_output_2d_on_each_pe )  THEN
              CALL check_open( 21 )
@@ -433,15 +427,24 @@
                 IF ( mode == 'xy' )  level_z = zu
 
-             CASE ( 'c_liq*_xy' )        ! 2d-array
-                IF ( av == 0 )  THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = c_liq(j,i) * c_veg(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = c_liq_av(j,i)
+             CASE ( 'lpt_xy', 'lpt_xz', 'lpt_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => pt
+                ELSE
+                   to_be_resorted => lpt_av
+                ENDIF
+                IF ( mode == 'xy' )  level_z = zu
+
+             CASE ( 'lwp*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = SUM( ql(nzb:nzt,j,i) *          &
+                                                    dzw(1:nzt+1) )
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = lwp_av(j,i)
                       ENDDO
                    ENDDO
@@ -450,131 +453,4 @@
                 two_d = .TRUE.
                 level_z(nzb+1) = zu(nzb+1)
-
-             CASE ( 'c_soil*_xy' )        ! 2d-array
-                IF ( av == 0 )  THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = 1.0_wp - c_veg(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = c_soil_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
-
-             CASE ( 'c_veg*_xy' )        ! 2d-array
-                IF ( av == 0 )  THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = c_veg(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = c_veg_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
-
-             CASE ( 'ghf_eb*_xy' )        ! 2d-array
-                IF ( av == 0 )  THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = ghf_eb(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = ghf_eb_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
-
-             CASE ( 'lai*_xy' )        ! 2d-array
-                IF ( av == 0 )  THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = lai(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = lai_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
-
-             CASE ( 'lpt_xy', 'lpt_xz', 'lpt_yz' )
-                IF ( av == 0 )  THEN
-                   to_be_resorted => pt
-                ELSE
-                   to_be_resorted => lpt_av
-                ENDIF
-                IF ( mode == 'xy' )  level_z = zu
-
-             CASE ( 'lwp*_xy' )        ! 2d-array
-                IF ( av == 0 )  THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = SUM( ql(nzb:nzt,j,i) *          &
-                                                    dzw(1:nzt+1) )
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = lwp_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
-
-             CASE ( 'm_liq_eb*_xy' )        ! 2d-array
-                IF ( av == 0 )  THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = m_liq_eb(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = m_liq_eb_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
-
-             CASE ( 'm_soil_xy', 'm_soil_xz', 'm_soil_yz' )
-                nzb_do = nzb_soil
-                nzt_do = nzt_soil
-                IF ( av == 0 )  THEN
-                   to_be_resorted => m_soil
-                ELSE
-                   to_be_resorted => m_soil_av
-                ENDIF
-                IF ( mode == 'xy' )  level_z = zs
 
              CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
@@ -859,9 +735,24 @@
                 level_z(nzb+1) = zu(nzb+1)
 
-             CASE ( 'qsws_eb*_xy' )        ! 2d-array
+             CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                   resorted = .TRUE.
+                ELSE
+                   to_be_resorted => qv_av
+                ENDIF
+                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) =  qsws_eb(j,i)
+                         local_pf(i,j,nzb+1) =  rad_net(j,i)
                       ENDDO
                    ENDDO
@@ -869,5 +760,5 @@
                    DO  i = nxlg, nxrg
                       DO  j = nysg, nyng 
-                         local_pf(i,j,nzb+1) =  qsws_eb_av(j,i)
+                         local_pf(i,j,nzb+1) =  rad_net_av(j,i)
                       ENDDO
                    ENDDO
@@ -877,91 +768,4 @@
                 level_z(nzb+1) = zu(nzb+1)
 
-             CASE ( 'qsws_liq_eb*_xy' )        ! 2d-array
-                IF ( av == 0 ) THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) =  qsws_liq_eb(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng 
-                         local_pf(i,j,nzb+1) =  qsws_liq_eb_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
-
-             CASE ( 'qsws_soil_eb*_xy' )        ! 2d-array
-                IF ( av == 0 ) THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) =  qsws_soil_eb(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng 
-                         local_pf(i,j,nzb+1) =  qsws_soil_eb_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
-
-             CASE ( 'qsws_veg_eb*_xy' )        ! 2d-array
-                IF ( av == 0 ) THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) =  qsws_veg_eb(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng 
-                         local_pf(i,j,nzb+1) =  qsws_veg_eb_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
-
-             CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
-                IF ( av == 0 )  THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         DO  k = nzb, nzt+1
-                            local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
-                         ENDDO
-                      ENDDO
-                   ENDDO
-                   resorted = .TRUE.
-                ELSE
-                   to_be_resorted => qv_av
-                ENDIF
-                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' )
@@ -1035,40 +839,4 @@
                    to_be_resorted => rho_av
                 ENDIF
-
-             CASE ( 'r_a*_xy' )        ! 2d-array
-                IF ( av == 0 )  THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = r_a(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = r_a_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
-
-             CASE ( 'r_s*_xy' )        ! 2d-array
-                IF ( av == 0 )  THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = r_s(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) = r_s_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
 
              CASE ( 's_xy', 's_xz', 's_yz' )
@@ -1097,22 +865,4 @@
                       DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) =  shf_av(j,i)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                resorted = .TRUE.
-                two_d = .TRUE.
-                level_z(nzb+1) = zu(nzb+1)
-
-             CASE ( 'shf_eb*_xy' )        ! 2d-array
-                IF ( av == 0 ) THEN
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) =  shf_eb(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) =  shf_eb_av(j,i)
                       ENDDO
                    ENDDO
@@ -1158,14 +908,4 @@
                 level_z(nzb+1) = zu(nzb+1)
 
-             CASE ( 't_soil_xy', 't_soil_xz', 't_soil_yz' )
-                nzb_do = nzb_soil
-                nzt_do = nzt_soil
-                IF ( av == 0 )  THEN
-                   to_be_resorted => t_soil
-                ELSE
-                   to_be_resorted => t_soil_av
-                ENDIF
-                IF ( mode == 'xy' )  level_z = zs
-
              CASE ( 'u_xy', 'u_xz', 'u_yz' )
                 IF ( av == 0 )  THEN
@@ -1285,8 +1025,19 @@
 
              CASE DEFAULT
+
+!
+!--             Land surface model quantity
+                IF ( land_surface )  THEN
+                   CALL lsm_data_output_2d( av, do2d(av,if), found, grid, mode,&
+                                            local_pf, two_d, nzb_do, nzt_do )
+                ENDIF
+
 !
 !--             User defined quantity
-                CALL user_data_output_2d( av, do2d(av,if), found, grid,        &
-                                          local_pf, two_d, nzb_do, nzt_do )
+                IF ( .NOT. found )  THEN
+                   CALL user_data_output_2d( av, do2d(av,if), found, grid,     &
+                                             local_pf, two_d, nzb_do, nzt_do )
+                ENDIF
+
                 resorted = .TRUE.
 
Index: /palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_3d.f90	(revision 1971)
+++ /palm/trunk/SOURCE/data_output_3d.f90	(revision 1972)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Output of land surface quantities is now done directly in the respective module.
+! Unnecessary directive __parallel removed.
 ! 
 ! Former revisions:
@@ -144,5 +145,5 @@
     
     USE land_surface_model_mod,                                                &
-        ONLY: m_soil, m_soil_av, nzb_soil, nzt_soil, t_soil, t_soil_av
+        ONLY: land_surface, lsm_data_output_3d, nzb_soil, nzt_soil
 
 #if defined( __netcdf )
@@ -272,18 +273,4 @@
              ELSE
                 to_be_resorted => lpt_av
-             ENDIF
-
-          CASE ( 'm_soil' )
-             nzb_do = nzb_soil
-             nzt_do = nzt_soil
-!
-!--          For soil model quantities, it is required to re-allocate local_pf
-             DEALLOCATE ( local_pf )
-             ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )
-
-             IF ( av == 0 )  THEN
-                to_be_resorted => m_soil
-             ELSE
-                to_be_resorted => m_soil_av
              ENDIF
 
@@ -581,18 +568,4 @@
              ENDIF
 
-          CASE ( 't_soil' )
-             nzb_do = nzb_soil
-             nzt_do = nzt_soil
-!
-!--          For soil model quantities, it is required to re-allocate local_pf
-             DEALLOCATE ( local_pf )
-             ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )
-
-             IF ( av == 0 )  THEN
-                to_be_resorted => t_soil
-             ELSE
-                to_be_resorted => t_soil_av
-             ENDIF
-
           CASE ( 'u' )
              IF ( av == 0 )  THEN
@@ -624,9 +597,27 @@
 
           CASE DEFAULT
+
+!
+!--          Land surface quantity
+             IF ( land_surface )  THEN
+!
+!--             For soil model quantities, it is required to re-allocate local_pf
+                nzb_do = nzb_soil
+                nzt_do = nzt_soil
+
+                DEALLOCATE ( local_pf )
+                ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )
+
+                CALL lsm_data_output_3d( av, do3d(av,if), found, local_pf )
+                resorted = .TRUE.
+             ENDIF
+
 !
 !--          User defined quantity
-             CALL user_data_output_3d( av, do3d(av,if), found, local_pf,       &
-                                       nzb_do, nzt_do )
-             resorted = .TRUE.
+             IF ( .NOT. found )  THEN
+                CALL user_data_output_3d( av, do3d(av,if), found, local_pf,    &
+                                          nzb_do, nzt_do )
+                resorted = .TRUE.
+             ENDIF
 
              IF ( .NOT. found )  THEN
@@ -667,7 +658,7 @@
                 WRITE ( 30 )  local_pf(:,:,nzb_do:nzt_do)
              ENDIF
-#if defined( __parallel )
+
              CALL MPI_BARRIER( comm2d, ierr )
-#endif
+
           ENDDO
 
Index: /palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 1971)
+++ /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 1972)
@@ -19,5 +19,7 @@
 ! Current revisions:
 ! -----------------
-! 
+! Further modularization: output of cross sections and 3D data is now done in this
+! module. Moreover, restart data is written and read directly within this module.
+!
 ! 
 ! Former revisions:
@@ -549,6 +551,8 @@
 !-- Public functions
     PUBLIC lsm_check_data_output, lsm_check_data_output_pr,                    &
-           lsm_check_parameters, lsm_energy_balance, lsm_header, lsm_init,     &
-           lsm_init_arrays, lsm_parin, lsm_soil_model, lsm_swap_timelevel 
+           lsm_check_parameters, lsm_define_netcdf_grid, lsm_3d_data_averaging,& 
+           lsm_data_output_2d, lsm_data_output_3d, lsm_energy_balance,         &
+           lsm_header, lsm_init, lsm_init_arrays, lsm_parin, lsm_soil_model,   &
+           lsm_swap_timelevel, lsm_read_restart_data, lsm_last_actions
 !
 !-- Public parameters, constants and initial values
@@ -561,12 +565,10 @@
 !
 !-- Public 2D output variables
-    PUBLIC c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb, ghf_eb_av,     &
-           lai, lai_av, qsws_eb, qsws_eb_av, qsws_liq_eb, qsws_liq_eb_av,      &
-           qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av,         &
-           r_a, r_a_av, r_s, r_s_av, shf_eb, shf_eb_av
+    PUBLIC ghf_eb, qsws_eb, qsws_liq_eb, qsws_soil_eb,qsws_veg_eb, r_a, r_s,   &
+           shf_eb
 
 !
 !-- Public prognostic variables
-    PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
+    PUBLIC m_soil, t_soil
 
 
@@ -583,4 +585,20 @@
     END INTERFACE lsm_check_parameters
     
+    INTERFACE lsm_3d_data_averaging
+       MODULE PROCEDURE lsm_3d_data_averaging
+    END INTERFACE lsm_3d_data_averaging
+
+    INTERFACE lsm_data_output_2d
+       MODULE PROCEDURE lsm_data_output_2d
+    END INTERFACE lsm_data_output_2d
+
+    INTERFACE lsm_data_output_3d
+       MODULE PROCEDURE lsm_data_output_3d
+    END INTERFACE lsm_data_output_3d
+
+    INTERFACE lsm_define_netcdf_grid
+       MODULE PROCEDURE lsm_define_netcdf_grid
+    END INTERFACE lsm_define_netcdf_grid
+
     INTERFACE lsm_energy_balance
        MODULE PROCEDURE lsm_energy_balance
@@ -611,4 +629,13 @@
     END INTERFACE lsm_swap_timelevel
 
+    INTERFACE lsm_read_restart_data
+       MODULE PROCEDURE lsm_read_restart_data
+    END INTERFACE lsm_read_restart_data
+
+
+    INTERFACE lsm_last_actions
+       MODULE PROCEDURE lsm_last_actions
+    END INTERFACE lsm_last_actions
+
  CONTAINS
 
@@ -627,5 +654,5 @@
 
        CHARACTER (LEN=*) ::  unit     !< 
-       CHARACTER (LEN=*) ::  var      !<
+       CHARACTER (LEN=*) ::  var !<
 
        INTEGER(iwp) :: i
@@ -2415,4 +2442,1093 @@
 
 
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Soubroutine for averaging 3D data
+!------------------------------------------------------------------------------!
+SUBROUTINE lsm_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 ( 'c_liq*' )
+                IF ( .NOT. ALLOCATED( c_liq_av ) )  THEN
+                   ALLOCATE( c_liq_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                c_liq_av = 0.0_wp
+
+             CASE ( 'c_soil*' )
+                IF ( .NOT. ALLOCATED( c_soil_av ) )  THEN
+                   ALLOCATE( c_soil_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                c_soil_av = 0.0_wp
+
+             CASE ( 'c_veg*' )
+                IF ( .NOT. ALLOCATED( c_veg_av ) )  THEN
+                   ALLOCATE( c_veg_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                c_veg_av = 0.0_wp
+
+             CASE ( 'ghf_eb*' )
+                IF ( .NOT. ALLOCATED( ghf_eb_av ) )  THEN
+                   ALLOCATE( ghf_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                ghf_eb_av = 0.0_wp
+
+             CASE ( 'lai*' )
+                IF ( .NOT. ALLOCATED( lai_av ) )  THEN
+                   ALLOCATE( lai_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                lai_av = 0.0_wp
+
+             CASE ( 'm_liq_eb*' )
+                IF ( .NOT. ALLOCATED( m_liq_eb_av ) )  THEN
+                   ALLOCATE( m_liq_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                m_liq_eb_av = 0.0_wp
+
+             CASE ( 'm_soil' )
+                IF ( .NOT. ALLOCATED( m_soil_av ) )  THEN
+                   ALLOCATE( m_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                m_soil_av = 0.0_wp
+
+             CASE ( 'qsws_eb*' )
+                IF ( .NOT. ALLOCATED( qsws_eb_av ) )  THEN
+                   ALLOCATE( qsws_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                qsws_eb_av = 0.0_wp
+
+             CASE ( 'qsws_liq_eb*' )
+                IF ( .NOT. ALLOCATED( qsws_liq_eb_av ) )  THEN
+                   ALLOCATE( qsws_liq_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                qsws_liq_eb_av = 0.0_wp
+
+             CASE ( 'qsws_soil_eb*' )
+                IF ( .NOT. ALLOCATED( qsws_soil_eb_av ) )  THEN
+                   ALLOCATE( qsws_soil_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                qsws_soil_eb_av = 0.0_wp
+
+             CASE ( 'qsws_veg_eb*' )
+                IF ( .NOT. ALLOCATED( qsws_veg_eb_av ) )  THEN
+                   ALLOCATE( qsws_veg_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                qsws_veg_eb_av = 0.0_wp
+
+             CASE ( 'r_a*' )
+                IF ( .NOT. ALLOCATED( r_a_av ) )  THEN
+                   ALLOCATE( r_a_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                r_a_av = 0.0_wp
+
+             CASE ( 'r_s*' )
+                IF ( .NOT. ALLOCATED( r_s_av ) )  THEN
+                   ALLOCATE( r_s_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                r_s_av = 0.0_wp
+
+             CASE ( 'shf_eb*' )
+                IF ( .NOT. ALLOCATED( shf_eb_av ) )  THEN
+                   ALLOCATE( shf_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                shf_eb_av = 0.0_wp
+
+             CASE ( 't_soil' )
+                IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN
+                   ALLOCATE( t_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                t_soil_av = 0.0_wp
+
+          CASE DEFAULT
+             CONTINUE
+
+       END SELECT
+
+    ELSEIF ( mode == 'sum' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+          CASE ( 'c_liq*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   c_liq_av(j,i) = c_liq_av(j,i) + c_liq(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'c_soil*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   c_soil_av(j,i) = c_soil_av(j,i) + (1.0 - c_veg(j,i))
+                ENDDO
+             ENDDO
+
+          CASE ( 'c_veg*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   c_veg_av(j,i) = c_veg_av(j,i) + c_veg(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'ghf_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   ghf_eb_av(j,i) = ghf_eb_av(j,i) + ghf_eb(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'lai*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   lai_av(j,i) = lai_av(j,i) + lai(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'm_liq_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   m_liq_eb_av(j,i) = m_liq_eb_av(j,i) + m_liq_eb(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'm_soil' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb_soil, nzt_soil
+                      m_soil_av(k,j,i) = m_soil_av(k,j,i) + m_soil(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'qsws_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_eb_av(j,i) = qsws_eb_av(j,i) + qsws_eb(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'qsws_liq_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) + qsws_liq_eb(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'qsws_soil_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) + qsws_soil_eb(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'qsws_veg_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) + qsws_veg_eb(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'r_a*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   r_a_av(j,i) = r_a_av(j,i) + r_a(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'r_s*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   r_s_av(j,i) = r_s_av(j,i) + r_s(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'shf_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   shf_eb_av(j,i) = shf_eb_av(j,i) + shf_eb(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 't_soil' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb_soil, nzt_soil
+                      t_soil_av(k,j,i) = t_soil_av(k,j,i) + t_soil(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE DEFAULT
+             CONTINUE
+
+       END SELECT
+
+    ELSEIF ( mode == 'average' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+          CASE ( 'c_liq*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   c_liq_av(j,i) = c_liq_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 'c_soil*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   c_soil_av(j,i) = c_soil_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 'c_veg*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   c_veg_av(j,i) = c_veg_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 'ghf_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   ghf_eb_av(j,i) = ghf_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+         CASE ( 'lai*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   lai_av(j,i) = lai_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 'm_liq_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   m_liq_eb_av(j,i) = m_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 'm_soil' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb_soil, nzt_soil
+                      m_soil_av(k,j,i) = m_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'qsws_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_eb_av(j,i) = qsws_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 'qsws_liq_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 'qsws_soil_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 'qsws_veg_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 'r_a*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   r_a_av(j,i) = r_a_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 'r_s*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   r_s_av(j,i) = r_s_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 't_soil' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb_soil, nzt_soil
+                      t_soil_av(k,j,i) = t_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+       END SELECT
+
+    ENDIF
+
+END SUBROUTINE lsm_3d_data_averaging
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Soubroutine defines 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 )
+    
+        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
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Soubroutine defines 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
+
+    CHARACTER (LEN=*) ::  grid     !< 
+    CHARACTER (LEN=*) ::  mode     !< 
+    CHARACTER (LEN=*) ::  variable !< 
+
+    INTEGER(iwp) ::  av !< 
+    INTEGER(iwp) ::  i  !< 
+    INTEGER(iwp) ::  j  !< 
+    INTEGER(iwp) ::  k  !< 
+    INTEGER(iwp) ::  nzb_do  !< 
+    INTEGER(iwp) ::  nzt_do  !< 
+
+    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 ( 'c_liq*_xy' )        ! 2d-array
+          IF ( av == 0 )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = c_liq(j,i) * c_veg(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = c_liq_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+       CASE ( 'c_soil*_xy' )        ! 2d-array
+          IF ( av == 0 )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = 1.0_wp - c_veg(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = c_soil_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+       CASE ( 'c_veg*_xy' )        ! 2d-array
+          IF ( av == 0 )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = c_veg(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = c_veg_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+       CASE ( 'ghf_eb*_xy' )        ! 2d-array
+          IF ( av == 0 )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = ghf_eb(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = ghf_eb_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+       CASE ( 'lai*_xy' )        ! 2d-array
+          IF ( av == 0 )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = lai(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = lai_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+       CASE ( 'm_liq_eb*_xy' )        ! 2d-array
+          IF ( av == 0 )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = m_liq_eb(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = m_liq_eb_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+       CASE ( 'm_soil_xy', 'm_soil_xz', 'm_soil_yz' )
+          IF ( av == 0 )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO k = nzb_soil, nzt_soil
+                      local_pf(i,j,k) = m_soil(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO k = nzb_soil, nzt_soil
+                      local_pf(i,j,k) = m_soil_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+          nzb_do = nzb_soil
+          nzt_do = nzt_soil
+
+          IF ( mode == 'xy' ) grid = 'zs'
+
+       CASE ( 'qsws_eb*_xy' )        ! 2d-array
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) =  qsws_eb(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   local_pf(i,j,nzb+1) =  qsws_eb_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+       CASE ( 'qsws_liq_eb*_xy' )        ! 2d-array
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) =  qsws_liq_eb(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   local_pf(i,j,nzb+1) =  qsws_liq_eb_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+       CASE ( 'qsws_soil_eb*_xy' )        ! 2d-array
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) =  qsws_soil_eb(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   local_pf(i,j,nzb+1) =  qsws_soil_eb_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+       CASE ( 'qsws_veg_eb*_xy' )        ! 2d-array
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) =  qsws_veg_eb(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   local_pf(i,j,nzb+1) =  qsws_veg_eb_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+
+       CASE ( 'r_a*_xy' )        ! 2d-array
+          IF ( av == 0 )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = r_a(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = r_a_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+       CASE ( 'r_s*_xy' )        ! 2d-array
+          IF ( av == 0 )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = r_s(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = r_s_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+       CASE ( 'shf_eb*_xy' )        ! 2d-array
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) =  shf_eb(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) =  shf_eb_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+       CASE ( 't_soil_xy', 't_soil_xz', 't_soil_yz' )
+          IF ( av == 0 )  THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO k = nzb_soil, nzt_soil
+                      local_pf(i,j,k) = t_soil(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO k = nzb_soil, nzt_soil
+                      local_pf(i,j,k) = t_soil_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+          nzb_do = nzb_soil
+          nzt_do = nzt_soil
+
+          IF ( mode == 'xy' )  grid = 'zs'
+
+       CASE DEFAULT
+          found = .FALSE.
+          grid  = 'none'
+
+    END SELECT
+ 
+ END SUBROUTINE lsm_data_output_2d
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Soubroutine defines 3D output variables
+!------------------------------------------------------------------------------!
+ SUBROUTINE lsm_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_soil:nzt_soil) ::  local_pf !< 
+
+
+    found = .TRUE.
+
+
+    SELECT CASE ( TRIM( variable ) )
+
+
+      CASE ( 'm_soil' )
+
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb_soil, nzt_soil
+                     local_pf(i,j,k) = m_soil(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb_soil, nzt_soil
+                     local_pf(i,j,k) = m_soil_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 't_soil' )
+
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb_soil, nzt_soil
+                     local_pf(i,j,k) = t_soil(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb_soil, nzt_soil
+                     local_pf(i,j,k) = t_soil_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+
+       CASE DEFAULT
+          found = .FALSE.
+
+    END SELECT
+
+
+ END SUBROUTINE lsm_data_output_3d
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Write restart data for land surface model
+!------------------------------------------------------------------------------!
+ SUBROUTINE lsm_last_actions
+ 
+
+    USE control_parameters
+        
+    USE kinds
+
+    IMPLICIT NONE
+
+    IF ( write_binary(1:4) == 'true' )  THEN
+       IF ( ALLOCATED( c_liq_av ) )  THEN
+          WRITE ( 14 )  'c_liq_av            ';  WRITE ( 14 ) c_liq_av
+       ENDIF
+       IF ( ALLOCATED( c_soil_av ) )  THEN
+          WRITE ( 14 )  'c_soil_av           ';  WRITE ( 14 ) c_soil_av
+       ENDIF
+       IF ( ALLOCATED( c_veg_av ) )  THEN
+          WRITE ( 14 )  'c_veg_av            ';  WRITE ( 14 ) c_veg_av
+       ENDIF
+       IF ( ALLOCATED( ghf_eb_av ) )  THEN
+          WRITE ( 14 )  'ghf_eb_av           ';  WRITE ( 14 )  ghf_eb_av
+       ENDIF
+       IF ( ALLOCATED( lai_av ) )  THEN
+          WRITE ( 14 )  'lai_av              ';  WRITE ( 14 )  lai_av
+       ENDIF
+       WRITE ( 14 )  'm_liq_eb            ';  WRITE ( 14 )  m_liq_eb
+       IF ( ALLOCATED( m_liq_eb_av ) )  THEN
+          WRITE ( 14 )  'm_liq_eb_av         ';  WRITE ( 14 )  m_liq_eb_av
+       ENDIF
+       WRITE ( 14 )  'm_soil              ';  WRITE ( 14 )  m_soil
+       IF ( ALLOCATED( m_soil_av ) )  THEN
+          WRITE ( 14 )  'm_soil_av           ';  WRITE ( 14 )  m_soil_av
+       ENDIF
+       IF ( ALLOCATED( qsws_eb_av ) )  THEN
+          WRITE ( 14 )  'qsws_eb_av          ';  WRITE ( 14 )  qsws_eb_av
+       ENDIF   
+       IF ( ALLOCATED( qsws_liq_eb_av ) )  THEN
+          WRITE ( 14 )  'qsws_liq_eb_av      ';  WRITE ( 14 )  qsws_liq_eb_av
+       ENDIF  
+       IF ( ALLOCATED( qsws_soil_eb_av ) )  THEN
+          WRITE ( 14 )  'qsws_soil_eb_av     ';  WRITE ( 14 )  qsws_soil_eb_av
+       ENDIF 
+       IF ( ALLOCATED( qsws_veg_eb_av ) )  THEN
+          WRITE ( 14 )  'qsws_veg_eb_av      ';  WRITE ( 14 )  qsws_veg_eb_av
+       ENDIF 
+       IF ( ALLOCATED( shf_eb_av ) )  THEN
+          WRITE ( 14 )  'shf_eb_av           ';  WRITE ( 14 )  shf_eb_av
+       ENDIF
+       WRITE ( 14 )  't_soil              ';  WRITE ( 14 )  t_soil
+       IF ( ALLOCATED( t_soil_av ) )  THEN
+          WRITE ( 14 )  't_soil_av           ';  WRITE ( 14 )  t_soil_av
+       ENDIF
+
+       WRITE ( 14 )  '*** end lsm ***     '
+
+    ENDIF
+
+ END SUBROUTINE lsm_last_actions
+
+
+SUBROUTINE 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,      &
+                                     tmp_2d )
+ 
+
+    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_soil:nzt_soil+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
+          tmp_3d   !< 
+
+    REAL(wp),                                                                  &
+       DIMENSION(nzb_soil:nzt_soil,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
+          tmp_3d2   !< 
+
+
+   IF ( initializing_actions == 'read_restart_data' )  THEN
+      READ ( 13 )  field_char
+
+      DO  WHILE ( TRIM( field_char ) /= '*** end lsm ***' )
+
+         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 ( 'c_liq_av' )
+                   IF ( .NOT. ALLOCATED( c_liq_av ) )  THEN
+                      ALLOCATE( c_liq_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   c_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'c_soil_av' )
+                   IF ( .NOT. ALLOCATED( c_soil_av ) )  THEN
+                      ALLOCATE( c_soil_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   c_soil_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'c_veg_av' )
+                   IF ( .NOT. ALLOCATED( c_veg_av ) )  THEN
+                      ALLOCATE( c_veg_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   c_veg_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'ghf_eb_av' )
+                   IF ( .NOT. ALLOCATED( ghf_eb_av ) )  THEN
+                      ALLOCATE( ghf_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   ghf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'm_liq_eb' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   m_liq_eb(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =        &
+                                 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'lai_av' )
+                   IF ( .NOT. ALLOCATED( lai_av ) )  THEN
+                      ALLOCATE( lai_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   lai_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'm_liq_eb_av' )
+                   IF ( .NOT. ALLOCATED( m_liq_eb_av ) )  THEN
+                      ALLOCATE( m_liq_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   m_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
+                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'm_soil' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d2(:,:,:)
+                   m_soil(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
+                          tmp_3d2(nzb_soil:nzt_soil,nysf-nbgp:nynf             &
+                          +nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'm_soil_av' )
+                   IF ( .NOT. ALLOCATED( m_soil_av ) )  THEN
+                      ALLOCATE( m_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d2(:,:,:)
+                   m_soil_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
+                                    tmp_3d2(nzb_soil:nzt_soil,nysf             &
+                                    -nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qsws_eb_av' )
+                   IF ( .NOT. ALLOCATED( qsws_eb_av ) )  THEN
+                      ALLOCATE( qsws_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qsws_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qsws_liq_eb_av' )
+                   IF ( .NOT. ALLOCATED( qsws_liq_eb_av ) )  THEN
+                      ALLOCATE( qsws_liq_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qsws_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                CASE ( 'qsws_soil_eb_av' )
+                   IF ( .NOT. ALLOCATED( qsws_soil_eb_av ) )  THEN
+                      ALLOCATE( qsws_soil_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qsws_soil_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qsws_veg_eb_av' )
+                   IF ( .NOT. ALLOCATED( qsws_veg_eb_av ) )  THEN
+                      ALLOCATE( qsws_veg_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qsws_veg_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =  &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'shf_eb_av' )
+                   IF ( .NOT. ALLOCATED( shf_eb_av ) )  THEN
+                      ALLOCATE( shf_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   shf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 't_soil' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   t_soil(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
+                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,               &
+                                                nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 't_soil_av' )
+                   IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN
+                      ALLOCATE( t_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d2(:,:,:)
+                   t_soil_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( 'lsm_read_restart_data', 'PA0441', 1, 2, 0, 6, &
+                                 0 )
+
+            END SELECT
+
+         ENDDO
+
+         READ ( 13 )  field_char
+
+      ENDDO
+   ENDIF
+
+ END SUBROUTINE lsm_read_restart_data
+
 !------------------------------------------------------------------------------!
 ! Description:
Index: /palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 1971)
+++ /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 1972)
@@ -19,5 +19,7 @@
 ! Current revisions:
 ! ------------------
-! 
+! Bugfix: wrong units for lsm quantities.
+! Definition of grids for land surface quantities is now done directly in the
+! respective module.
 ! 
 ! Former revisions:
@@ -208,7 +210,7 @@
              'm2/s2        ', 'm2/s2        ', 'kg m/s       ',                &
              'm            ', 'kg/kg        ', 'kg m/(kg s)  ',                &
-             'kg/kg        ', '             ', '             ',                &
-             '             ', '             ', '             ',                &
-             'W/m2         ', 's/m          ', '             ',                &
+             'kg/kg        ', 'W/m2         ', 'W/m2         ',                &
+             'W/m2         ', 'W/m2         ', 'W/m2         ',                &
+             'W/m2         ', 's/m          ', 's/m          ',                &
              'W/m2         ', 'W/m2         ', 'W/m2         ',                &
              'W/m2         ', 'W/m2         ', '             ',                &
@@ -400,5 +402,5 @@
 
     USE land_surface_model_mod,                                                &
-        ONLY: land_surface, nzb_soil, nzt_soil, nzs, zs
+        ONLY: land_surface, lsm_define_netcdf_grid, nzb_soil, nzt_soil, nzs, zs
 
     USE pegrid
@@ -758,17 +760,21 @@
                    grid_y = 'y'
                    grid_z = 'zw'
-!
-!--             soil grid
-                CASE ( 'm_soil', 't_soil' )
-
-                   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( domask(mid,av,i), found,    &
+                                                   grid_x, grid_y, grid_z )
+                   ENDIF
+
 !
 !--                Check for user-defined quantities
-                   CALL user_define_netcdf_grid( domask(mid,av,i), found,      &
-                                                 grid_x, grid_y, grid_z )
+                   IF ( .NOT. found )  THEN
+                      CALL user_define_netcdf_grid( domask(mid,av,i), found,   &
+                                                    grid_x, grid_y, grid_z )
+                   ENDIF
 
                    IF ( .NOT. found )  THEN
@@ -1258,17 +1264,20 @@
                    grid_y = 'y'
                    grid_z = 'zw'
-!
-!--             soil grid
-                CASE ( 'm_soil', 't_soil' )
-
-                   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( do3d(av,i), found, grid_x,  &
+                                                    grid_y, grid_z )
+                   ENDIF
+                   
 !--                Check for user-defined quantities
-                   CALL user_define_netcdf_grid( do3d(av,i), found, grid_x, &
-                                                 grid_y, grid_z )
+                   IF ( .NOT. found )  THEN
+                      CALL user_define_netcdf_grid( do3d(av,i), found, grid_x, &
+                                                    grid_y, grid_z )
+                   ENDIF
 
                    IF ( .NOT. found )  THEN
@@ -1992,5 +2001,5 @@
 
 !
-!--             Write zs data
+!--          Write zs data
              IF ( land_surface )  THEN
                 ns_do = 0
Index: /palm/trunk/SOURCE/palm.f90
===================================================================
--- /palm/trunk/SOURCE/palm.f90	(revision 1971)
+++ /palm/trunk/SOURCE/palm.f90	(revision 1972)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added call to lsm_last_actions for binary output of land surface model data
 ! 
 ! Former revisions:
@@ -161,4 +161,7 @@
 
     USE kinds
+
+    USE land_surface_model_mod,                                                &
+        ONLY:  land_surface, lsm_last_actions
 
     USE ls_forcing_mod,                                                        &
@@ -442,5 +445,5 @@
     IF ( myid == 0 )  CALL header
 !
-!-- If required, final user-defined actions, and
+!-- If required, final land surface and user-defined actions, and
 !-- last actions on the open files and close files. Unit 14 was opened
 !-- in write_3d_binary but it is closed here, to allow writing on this
@@ -449,4 +452,7 @@
     DO  i = 0, io_blocks-1
        IF ( i == io_group )  THEN
+          IF ( land_surface )  THEN
+             CALL lsm_last_actions
+          ENDIF
           CALL user_last_actions
           IF ( write_binary(1:4) == 'true' )  CALL close_file( 14 )
Index: /palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/read_3d_binary.f90	(revision 1971)
+++ /palm/trunk/SOURCE/read_3d_binary.f90	(revision 1972)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Land surface parts are now done in the respective module
 ! 
 ! Former revisions:
@@ -127,8 +127,5 @@
 
     USE land_surface_model_mod,                                                &
-        ONLY:  c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb_av, lai_av,        &
-               qsws_eb_av, qsws_liq_eb_av, qsws_soil_eb_av, qsws_veg_eb_av,    & 
-               land_surface, m_liq_eb, m_liq_eb_av, m_soil, m_soil_av,         &
-               nzb_soil, nzt_soil, shf_eb_av, t_soil, t_soil_av
+        ONLY:  land_surface, lsm_read_restart_data
 
     USE particle_attributes,                                                   &
@@ -204,6 +201,4 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  tmp_2d      !< temporary array for storing 2D data
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3d      !< temporary array for storing 3D data
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3d_soil1!< temporary array for storing 3D soil model data
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3d_soil2!< temporary array for storing 3D soil model data
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwul   !<
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwun   !<
@@ -330,5 +325,5 @@
 !--    First compare the version numbers
        READ ( 13 )  version_on_file
-       binary_version = '4.3'
+       binary_version = '4.4'
        IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
           WRITE( message_string, * ) 'version mismatch concerning data ',      &
@@ -405,15 +400,4 @@
                         nxl_on_file-nbgp:nxr_on_file+nbgp) )
 
-!
-!--    If the land surface model is used, allocate array for storing soil model
-!--    data
-       IF ( land_surface )  THEN
-          ALLOCATE( tmp_3d_soil1(nzb_soil:nzt_soil+1,                          &
-                                nys_on_file-nbgp:nyn_on_file+nbgp,             &
-                                nxl_on_file-nbgp:nxr_on_file+nbgp) )
-          ALLOCATE( tmp_3d_soil2(nzb_soil:nzt_soil,                            &
-                                nys_on_file-nbgp:nyn_on_file+nbgp,             &
-                                nxl_on_file-nbgp:nxr_on_file+nbgp) )
-       ENDIF
 !
 !--    Read arrays
@@ -444,28 +428,4 @@
              SELECT CASE ( TRIM( field_chr ) )
 
-                CASE ( 'c_liq_av' )
-                   IF ( .NOT. ALLOCATED( c_liq_av ) )  THEN
-                      ALLOCATE( c_liq_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   c_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
-                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'c_soil_av' )
-                   IF ( .NOT. ALLOCATED( c_soil_av ) )  THEN
-                      ALLOCATE( c_soil_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   c_soil_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
-                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'c_veg_av' )
-                   IF ( .NOT. ALLOCATED( c_veg_av ) )  THEN
-                      ALLOCATE( c_veg_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   c_veg_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
-                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
                 CASE ( 'e' )
                    IF ( k == 1 )  READ ( 13 )  tmp_3d
@@ -481,24 +441,8 @@
                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
-                CASE ( 'ghf_eb_av' )
-                   IF ( .NOT. ALLOCATED( ghf_eb_av ) )  THEN
-                      ALLOCATE( ghf_eb_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   ghf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
-                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
                 CASE ( 'iran' ) ! matching random numbers is still unresolved
                                 ! issue
                    IF ( k == 1 )  READ ( 13 )  iran, iran_part
 
-                CASE ( 'lai_av' )
-                   IF ( .NOT. ALLOCATED( lai_av ) )  THEN
-                      ALLOCATE( lai_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   lai_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
-                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
                 CASE ( 'kh' )
                    IF ( k == 1 )  READ ( 13 )  tmp_3d
@@ -526,31 +470,4 @@
                    lwp_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
                                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'm_liq_eb' )
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   m_liq_eb(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =        &
-                                 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-                CASE ( 'm_liq_eb_av' )
-                   IF ( .NOT. ALLOCATED( m_liq_eb_av ) )  THEN
-                      ALLOCATE( m_liq_eb_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   m_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
-                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'm_soil' )
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d_soil2(:,:,:)
-                   m_soil(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
-                          tmp_3d_soil2(nzb_soil:nzt_soil,nysf-nbgp:nynf        &
-                          +nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'm_soil_av' )
-                   IF ( .NOT. ALLOCATED( m_soil_av ) )  THEN
-                      ALLOCATE( m_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d_soil2(:,:,:)
-                   m_soil_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
-                                    tmp_3d_soil2(nzb_soil:nzt_soil,nysf        &
-                                    -nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'nr' )
@@ -769,35 +686,4 @@
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
                    qswst(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
-                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'qsws_eb_av' )
-                   IF ( .NOT. ALLOCATED( qsws_eb_av ) )  THEN
-                      ALLOCATE( qsws_eb_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF  
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   qsws_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
-                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'qsws_liq_eb_av' )
-                   IF ( .NOT. ALLOCATED( qsws_liq_eb_av ) )  THEN
-                      ALLOCATE( qsws_liq_eb_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF  
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   qsws_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
-                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-                CASE ( 'qsws_soil_eb_av' )
-                   IF ( .NOT. ALLOCATED( qsws_soil_eb_av ) )  THEN
-                      ALLOCATE( qsws_soil_eb_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF  
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   qsws_soil_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
-                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'qsws_veg_eb_av' )
-                   IF ( .NOT. ALLOCATED( qsws_veg_eb_av ) )  THEN
-                      ALLOCATE( qsws_veg_eb_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF  
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   qsws_veg_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -1033,12 +919,4 @@
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
                    shf_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
-                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'shf_eb_av' )
-                   IF ( .NOT. ALLOCATED( shf_eb_av ) )  THEN
-                      ALLOCATE( shf_eb_av(nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   shf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -1113,19 +991,4 @@
                    tswst(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 't_soil' )
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d_soil1
-                   t_soil(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
-                                   tmp_3d_soil1(:,nysf-nbgp:nynf+nbgp,         &
-                                                nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 't_soil_av' )
-                   IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN
-                      ALLOCATE( t_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-                   ENDIF
-                   IF ( k == 1 )  READ ( 13 )  tmp_3d_soil2(:,:,:)
-                   t_soil_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
-                                    tmp_3d_soil2(:,nysf-nbgp:nynf+nbgp,        &
-                                    nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'u' )
@@ -1407,4 +1270,11 @@
 
 !
+!--    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 )
+
+!
 !--    Read user-defined restart data
        CALL user_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file, &
@@ -1418,6 +1288,4 @@
 
        DEALLOCATE( tmp_2d, tmp_3d )
-       IF ( ALLOCATED( tmp_3d_soil1 ) )  DEALLOCATE( tmp_3d_soil1 )
-       IF ( ALLOCATED( tmp_3d_soil2 ) )  DEALLOCATE( tmp_3d_soil2 )
 
     ENDDO  ! loop over restart files
Index: /palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- /palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 1971)
+++ /palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 1972)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Land surface actions are now done directly in the respective module
 ! 
 ! Former revisions:
@@ -129,10 +129,5 @@
 
     USE land_surface_model_mod,                                                &
-        ONLY:  c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb,            &
-               ghf_eb_av, lai, lai_av, m_liq_eb, m_liq_eb_av, m_soil,          &
-               m_soil_av, nzb_soil, nzt_soil, qsws_eb, qsws_eb_av,             &
-               qsws_liq_eb, qsws_liq_eb_av, qsws_soil_eb, qsws_soil_eb_av,     &
-               qsws_veg_eb, qsws_veg_eb_av, shf_eb, shf_eb_av, r_a, r_a_av,    &
-               r_s, r_s_av, t_soil, t_soil_av
+        ONLY:  land_surface, lsm_3d_data_averaging
 
     USE particle_attributes,                                                   &
@@ -173,22 +168,4 @@
           SELECT CASE ( TRIM( doav(ii) ) )
 
-             CASE ( 'c_liq*' )
-                IF ( .NOT. ALLOCATED( c_liq_av ) )  THEN
-                   ALLOCATE( c_liq_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                c_liq_av = 0.0_wp
-
-             CASE ( 'c_soil*' )
-                IF ( .NOT. ALLOCATED( c_soil_av ) )  THEN
-                   ALLOCATE( c_soil_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                c_soil_av = 0.0_wp
-
-             CASE ( 'c_veg*' )
-                IF ( .NOT. ALLOCATED( c_veg_av ) )  THEN
-                   ALLOCATE( c_veg_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                c_veg_av = 0.0_wp
-
              CASE ( 'e' )
                 IF ( .NOT. ALLOCATED( e_av ) )  THEN
@@ -197,16 +174,4 @@
                 e_av = 0.0_wp
 
-             CASE ( 'ghf_eb*' )
-                IF ( .NOT. ALLOCATED( ghf_eb_av ) )  THEN
-                   ALLOCATE( ghf_eb_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                ghf_eb_av = 0.0_wp
-
-             CASE ( 'lai*' )
-                IF ( .NOT. ALLOCATED( lai_av ) )  THEN
-                   ALLOCATE( lai_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                lai_av = 0.0_wp
-
              CASE ( 'lpt' )
                 IF ( .NOT. ALLOCATED( lpt_av ) )  THEN
@@ -221,16 +186,4 @@
                 lwp_av = 0.0_wp
 
-             CASE ( 'm_liq_eb*' )
-                IF ( .NOT. ALLOCATED( m_liq_eb_av ) )  THEN
-                   ALLOCATE( m_liq_eb_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                m_liq_eb_av = 0.0_wp
-
-             CASE ( 'm_soil' )
-                IF ( .NOT. ALLOCATED( m_soil_av ) )  THEN
-                   ALLOCATE( m_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                m_soil_av = 0.0_wp
-
              CASE ( 'nr' )
                 IF ( .NOT. ALLOCATED( nr_av ) )  THEN
@@ -329,28 +282,4 @@
                 qsws_av = 0.0_wp
 
-             CASE ( 'qsws_eb*' )
-                IF ( .NOT. ALLOCATED( qsws_eb_av ) )  THEN
-                   ALLOCATE( qsws_eb_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                qsws_eb_av = 0.0_wp
-
-             CASE ( 'qsws_liq_eb*' )
-                IF ( .NOT. ALLOCATED( qsws_liq_eb_av ) )  THEN
-                   ALLOCATE( qsws_liq_eb_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                qsws_liq_eb_av = 0.0_wp
-
-             CASE ( 'qsws_soil_eb*' )
-                IF ( .NOT. ALLOCATED( qsws_soil_eb_av ) )  THEN
-                   ALLOCATE( qsws_soil_eb_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                qsws_soil_eb_av = 0.0_wp
-
-             CASE ( 'qsws_veg_eb*' )
-                IF ( .NOT. ALLOCATED( qsws_veg_eb_av ) )  THEN
-                   ALLOCATE( qsws_veg_eb_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                qsws_veg_eb_av = 0.0_wp
-
              CASE ( 'qv' )
                 IF ( .NOT. ALLOCATED( qv_av ) )  THEN
@@ -419,16 +348,4 @@
                 rho_av = 0.0_wp
 
-             CASE ( 'r_a*' )
-                IF ( .NOT. ALLOCATED( r_a_av ) )  THEN
-                   ALLOCATE( r_a_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                r_a_av = 0.0_wp
-
-             CASE ( 'r_s*' )
-                IF ( .NOT. ALLOCATED( r_s_av ) )  THEN
-                   ALLOCATE( r_s_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                r_s_av = 0.0_wp
-
              CASE ( 's' )
                 IF ( .NOT. ALLOCATED( s_av ) )  THEN
@@ -449,16 +366,4 @@
                 shf_av = 0.0_wp
 
-             CASE ( 'shf_eb*' )
-                IF ( .NOT. ALLOCATED( shf_eb_av ) )  THEN
-                   ALLOCATE( shf_eb_av(nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                shf_eb_av = 0.0_wp
-
-             CASE ( 't_soil' )
-                IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN
-                   ALLOCATE( t_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-                ENDIF
-                t_soil_av = 0.0_wp
-
              CASE ( 't*' )
                 IF ( .NOT. ALLOCATED( ts_av ) )  THEN
@@ -516,4 +421,11 @@
 
              CASE DEFAULT
+
+!
+!--             Land surface quantity
+                IF ( land_surface )  THEN
+                   CALL lsm_3d_data_averaging( 'allocate', doav(ii) )
+                ENDIF
+
 !
 !--             User-defined quantity
@@ -534,25 +446,4 @@
        SELECT CASE ( TRIM( doav(ii) ) )
 
-          CASE ( 'c_liq*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   c_liq_av(j,i) = c_liq_av(j,i) + c_liq(j,i)
-                ENDDO
-             ENDDO
-
-          CASE ( 'c_soil*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   c_soil_av(j,i) = c_soil_av(j,i) + (1.0 - c_veg(j,i))
-                ENDDO
-             ENDDO
-
-          CASE ( 'c_veg*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   c_veg_av(j,i) = c_veg_av(j,i) + c_veg(j,i)
-                ENDDO
-             ENDDO
-
           CASE ( 'e' )
              DO  i = nxlg, nxrg
@@ -561,18 +452,4 @@
                       e_av(k,j,i) = e_av(k,j,i) + e(k,j,i)
                    ENDDO
-                ENDDO
-             ENDDO
-
-          CASE ( 'ghf_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   ghf_eb_av(j,i) = ghf_eb_av(j,i) + ghf_eb(j,i)
-                ENDDO
-             ENDDO
-
-          CASE ( 'lai*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   lai_av(j,i) = lai_av(j,i) + lai(j,i)
                 ENDDO
              ENDDO
@@ -592,20 +469,4 @@
                    lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i)            &
                                                * dzw(1:nzt+1) ) * rho_surface
-                ENDDO
-             ENDDO
-
-          CASE ( 'm_liq_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   m_liq_eb_av(j,i) = m_liq_eb_av(j,i) + m_liq_eb(j,i)
-                ENDDO
-             ENDDO
-
-          CASE ( 'm_soil' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb_soil, nzt_soil
-                      m_soil_av(k,j,i) = m_soil_av(k,j,i) + m_soil(k,j,i)
-                   ENDDO
                 ENDDO
              ENDDO
@@ -782,32 +643,4 @@
              ENDDO
 
-          CASE ( 'qsws_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   qsws_eb_av(j,i) = qsws_eb_av(j,i) + qsws_eb(j,i)
-                ENDDO
-             ENDDO
-
-          CASE ( 'qsws_liq_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) + qsws_liq_eb(j,i)
-                ENDDO
-             ENDDO
-
-          CASE ( 'qsws_soil_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) + qsws_soil_eb(j,i)
-                ENDDO
-             ENDDO
-
-          CASE ( 'qsws_veg_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) + qsws_veg_eb(j,i)
-                ENDDO
-             ENDDO
-
           CASE ( 'qv' )
              DO  i = nxlg, nxrg
@@ -898,18 +731,4 @@
              ENDDO
 
-          CASE ( 'r_a*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   r_a_av(j,i) = r_a_av(j,i) + r_a(j,i)
-                ENDDO
-             ENDDO
-
-          CASE ( 'r_s*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   r_s_av(j,i) = r_s_av(j,i) + r_s(j,i)
-                ENDDO
-             ENDDO
-
           CASE ( 'rho' )
              DO  i = nxlg, nxrg
@@ -946,11 +765,4 @@
              ENDDO
 
-          CASE ( 'shf_eb*' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   shf_eb_av(j,i) = shf_eb_av(j,i) + shf_eb(j,i)
-                ENDDO
-             ENDDO
-
           CASE ( 'ssws*' )
              DO  i = nxlg, nxrg
@@ -967,13 +779,4 @@
              ENDDO
 
-          CASE ( 't_soil' )
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   DO  k = nzb_soil, nzt_soil
-                      t_soil_av(k,j,i) = t_soil_av(k,j,i) + t_soil(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
-
           CASE ( 'u' )
              DO  i = nxlg, nxrg
@@ -1041,4 +844,10 @@
 
           CASE DEFAULT
+!
+!--          Land surface quantity
+             IF ( land_surface )  THEN
+                CALL lsm_3d_data_averaging( 'sum', doav(ii) )
+             ENDIF
+
 !
 !--          User-defined quantity
Index: /palm/trunk/SOURCE/user_init_land_surface.f90
===================================================================
--- /palm/trunk/SOURCE/user_init_land_surface.f90	(revision 1971)
+++ /palm/trunk/SOURCE/user_init_land_surface.f90	(revision 1972)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Update of use statements
 ! 
 ! Former revisions:
@@ -41,6 +41,4 @@
  
 
-    USE arrays_3d
-    
     USE control_parameters
     
@@ -48,6 +46,11 @@
     
     USE kinds
+    
+    USE land_surface_model_mod
 
-    USE land_surface_model_mod
+    USE netcdf_interface,                                                      &
+        ONLY: dots_label, dots_unit, dots_num
+    
+    USE pegrid
     
     USE user
Index: /palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/write_3d_binary.f90	(revision 1971)
+++ /palm/trunk/SOURCE/write_3d_binary.f90	(revision 1972)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Land surface actions are now done directly in the respective module
 ! 
 ! Former revisions:
@@ -113,10 +113,4 @@
     USE indices,                                                               &
         ONLY:  nxl, nxr, nys, nyn, nzb, nzt
-
-    USE land_surface_model_mod,                                                &
-        ONLY:  c_liq_av, c_soil_av, c_veg_av, ghf_eb_av, lai_av,               &
-               qsws_eb_av, qsws_liq_eb_av, qsws_soil_eb_av, qsws_veg_eb_av,    & 
-               land_surface, m_liq_eb, m_liq_eb_av, m_soil, m_soil_av,         &
-               nzb_soil, nzt_soil, shf_eb_av, t_soil, t_soil_av
         
     USE particle_attributes,                                                   &
@@ -152,5 +146,5 @@
 !
 !-- Write arrays.
-    binary_version = '4.3'
+    binary_version = '4.4'
 
     WRITE ( 14 )  binary_version
@@ -164,16 +158,4 @@
 !--            adjusted accordingly.
 
-
-    IF ( land_surface )  THEN
-       IF ( ALLOCATED( c_liq_av ) )  THEN
-          WRITE ( 14 )  'c_liq_av                 ';  WRITE ( 14 ) c_liq_av
-       ENDIF
-       IF ( ALLOCATED( c_soil_av ) )  THEN
-          WRITE ( 14 )  'c_soil_av                ';  WRITE ( 14 ) c_soil_av
-       ENDIF
-       IF ( ALLOCATED( c_veg_av ) )  THEN
-          WRITE ( 14 )  'c_veg_av                 ';  WRITE ( 14 ) c_veg_av
-       ENDIF
-    ENDIF
     WRITE ( 14 )  'e                   ';  WRITE ( 14 )  e
     IF ( ALLOCATED( e_av ) )  THEN
@@ -181,16 +163,6 @@
     ENDIF
     WRITE ( 14 )  'iran                ';  WRITE ( 14 )  iran, iran_part
-    IF ( land_surface )  THEN
-       IF ( ALLOCATED( ghf_eb_av ) )  THEN
-          WRITE ( 14 )  'ghf_eb_av                ';  WRITE ( 14 )  ghf_eb_av
-       ENDIF
-    ENDIF
     WRITE ( 14 )  'kh                  ';  WRITE ( 14 )  kh
     WRITE ( 14 )  'km                  ';  WRITE ( 14 )  km
-    IF ( land_surface )  THEN
-       IF ( ALLOCATED( lai_av ) )  THEN
-          WRITE ( 14 )  'lai_av                   ';  WRITE ( 14 )  lai_av
-       ENDIF
-    ENDIF
     IF ( ALLOCATED( lpt_av ) )  THEN
        WRITE ( 14 )  'lpt_av              ';  WRITE ( 14 )  lpt_av
@@ -198,14 +170,4 @@
     IF ( ALLOCATED( lwp_av ) )  THEN
        WRITE ( 14 )  'lwp_av              ';  WRITE ( 14 )  lwp_av
-    ENDIF
-    IF ( land_surface )  THEN
-       WRITE ( 14 )  'm_liq_eb            ';  WRITE ( 14 )  m_liq_eb
-       IF ( ALLOCATED( m_liq_eb_av ) )  THEN
-          WRITE ( 14 )  'm_liq_eb_av              ';  WRITE ( 14 )  m_liq_eb_av
-       ENDIF
-       WRITE ( 14 )  'm_soil              ';  WRITE ( 14 )  m_soil
-       IF ( ALLOCATED( m_soil_av ) )  THEN
-          WRITE ( 14 )  'm_soil_av                ';  WRITE ( 14 )  m_soil_av
-       ENDIF
     ENDIF
     WRITE ( 14 )  'ol                  ';  WRITE ( 14 )  ol
@@ -289,18 +251,4 @@
        WRITE ( 14 )  'sswst               ';  WRITE ( 14 ) sswst
     ENDIF    
-    IF ( land_surface )  THEN
-       IF ( ALLOCATED( qsws_eb_av ) )  THEN
-          WRITE ( 14 )  'qsws_eb_av          ';  WRITE ( 14 )  qsws_eb_av
-       ENDIF   
-       IF ( ALLOCATED( qsws_liq_eb_av ) )  THEN
-          WRITE ( 14 )  'qsws_liq_eb_av      ';  WRITE ( 14 )  qsws_liq_eb_av
-       ENDIF  
-       IF ( ALLOCATED( qsws_soil_eb_av ) )  THEN
-          WRITE ( 14 )  'qsws_soil_eb_av     ';  WRITE ( 14 )  qsws_soil_eb_av
-       ENDIF 
-       IF ( ALLOCATED( qsws_veg_eb_av ) )  THEN
-          WRITE ( 14 )  'qsws_veg_eb_av      ';  WRITE ( 14 )  qsws_veg_eb_av
-       ENDIF 
-    ENDIF
     IF ( ALLOCATED( rad_net ) )  THEN
        WRITE ( 14 )  'rad_net             ';  WRITE ( 14 )  rad_net  
@@ -374,10 +322,4 @@
        WRITE ( 14 )  'saswst              ';  WRITE ( 14 )  saswst
     ENDIF
-    IF ( land_surface )  THEN
-       WRITE ( 14 )  't_soil              ';  WRITE ( 14 )  t_soil
-       IF ( ALLOCATED( t_soil_av ) )  THEN
-          WRITE ( 14 )  't_soil_av                ';  WRITE ( 14 )  t_soil_av
-       ENDIF
-    ENDIF
     IF ( ALLOCATED( ql_c_av ) )  THEN
        WRITE ( 14 )  'ql_c_av             ';  WRITE ( 14 )  ql_c_av
@@ -408,9 +350,4 @@
        WRITE ( 14 )  'shf_av              ';  WRITE ( 14 )  shf_av
     ENDIF
-    IF ( land_surface )  THEN
-       IF ( ALLOCATED( shf_eb_av ) )  THEN
-          WRITE ( 14 )  'shf_eb_av           ';  WRITE ( 14 )  shf_eb_av
-       ENDIF
-    ENDIF
     IF ( ALLOCATED( spectrum_x ) )  THEN
        WRITE ( 14 )  'spectrum_x          ';  WRITE ( 14 )  spectrum_x
