Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 1967)
+++ 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:
