Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 4256)
+++ palm/trunk/SOURCE/Makefile	(revision 4258)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Add dependency of land-surface model on pmc_handle_communicator and cpu_log
+# 
+# 4245 2019-09-30 08:40:37Z pavelkrc
 # Remove no longer needed dependencies on surface_mod
 # 
@@ -689,7 +692,9 @@
 	bulk_cloud_model_mod.o \
 	calc_mean_profile.o \
-	mod_kinds.o \
-	modules.o \
-	netcdf_data_input_mod.o \
+	cpulog_mod.o \
+	mod_kinds.o \
+	modules.o \
+	netcdf_data_input_mod.o \
+	pmc_handle_communicator_mod.o \
 	pmc_interface_mod.o \
 	radiation_model_mod.o \
Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 4256)
+++ palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 4258)
@@ -25,4 +25,14 @@
 ! -----------------
 ! $Id$
+! - Revise limitation for soil moisture in case it exceeds its saturation 
+!   value (J. Resler)
+! - Revise initialization of soil moisture and temperature in a nested run in 
+!   case dynamic input information is available. This case, the soil within 
+!   the child domains can be initialized separately. (J. Resler, M. Suehring)
+! - As part of this revision, migrate the netcdf input of soil temperature / 
+!   moisture to this module, as well as the routine to inter/extrapolate soil 
+!   profiles between different grids. 
+! 
+! 4251 2019-10-02 12:07:38Z maronga
 ! Bugfix: albedo_types for vegetation_type look-up table corrected.
 ! 
@@ -174,5 +184,7 @@
 
     USE control_parameters,                                                    &
-        ONLY:  cloud_droplets, coupling_start_time,                            &
+        ONLY:  cloud_droplets,                                                 &
+               coupling_char,                                                  &
+               coupling_start_time,                                            &
                debug_output, debug_output_timestep, debug_string,              &
                dt_3d,                                                          &
@@ -183,4 +195,7 @@
                surface_pressure, timestep_scheme, tsc,                         &
                time_since_reference_point
+               
+    USE cpulog,                                                                &
+        ONLY:  cpu_log, log_point_s
 
     USE indices,                                                               &
@@ -192,10 +207,20 @@
     USE netcdf_data_input_mod,                                                 &
         ONLY :  building_type_f,                                               &
+                char_fill,                                                     &
+                char_lod,                                                      &
+                check_existence,                                               &
+                close_input_file,                                              &
+                get_attribute,                                                 &
+                get_dimension_length,                                          &
+                get_variable,                                                  &
                 init_3d,                                                       &
                 input_file_dynamic,                                            &
                 input_pids_dynamic,                                            &
                 input_pids_static,                                             &
-                netcdf_data_input_interpolate,                                 &
-                netcdf_data_input_init_lsm,                                    &
+                inquire_num_variables,                                         &
+                inquire_variable_names,                                        &
+                num_var_pids,                                                  &
+                open_read_file,                                                &
+                pids_id,                                                       &
                 pavement_pars_f,                                               &
                 pavement_subsurface_pars_f,                                    &
@@ -205,8 +230,9 @@
                 soil_type_f,                                                   &
                 surface_fraction_f,                                            &
+                vars_pids,                                                     &
                 vegetation_pars_f,                                             &
                 vegetation_type_f,                                             &
                 water_pars_f,                                                  &
-                water_type_f
+                water_type_f             
 
     USE kinds
@@ -2276,11 +2302,11 @@
            ONLY:  nx, ny, topo_min_level
 
+       USE pmc_handle_communicator,                                            &
+        ONLY:  pmc_is_rootmodel
+           
        USE pmc_interface,                                                      &
            ONLY:  nested_run
     
        IMPLICIT NONE
-
-       LOGICAL      ::  init_msoil_from_parent   !< flag controlling initialization of soil moisture in nested child domains
-       LOGICAL      ::  init_tsoil_from_parent   !< flag controlling initialization of soil temperature in nested child domains
 
        INTEGER(iwp) ::  i                       !< running index
@@ -2295,7 +2321,17 @@
        INTEGER(iwp) ::  st                      !< soil-type index
        INTEGER(iwp) ::  n_soil_layers_total     !< temperature variable, stores the total number of soil layers + 4 
-
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  bound, bound_root_fr  !< temporary arrays for storing index bounds
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_soil_init !< temporary array used for averaging soil profiles
+#if defined( __parallel )
+       INTEGER(iwp) ::  nzs_root                !< number of soil layers in root domain (used in case soil data needs to be
+                                                !< transferred from root model to child models)   
+                                                
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  m_soil_root    !< domain-averaged soil moisture profile in root domain 
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_soil_root    !< domain-averaged soil temperature profile in root domain
+#endif
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  bound          !< temporary arrays for storing index bounds
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  bound_root_fr  !< temporary arrays for storing index bounds
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_soil_init   !< temporary array used for averaging soil profiles
+#if defined( __parallel )
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z_soil_root    !< vertical dimension of soil grid in root domain
+#endif
 
        IF ( debug_output )  CALL debug_message( 'lsm_init', 'start' )
@@ -4018,27 +4054,4 @@
        IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
 !
-!--       Read soil properties from dynamic input file. 
-          IF ( INDEX( initializing_actions, 'inifor' ) /= 0 )                  &
-             CALL netcdf_data_input_init_lsm
-!
-!--       In case no dynamic input is available for a child domain but root 
-!--       domain is initialized with dynamic input file, the different soil 
-!--       properties can lead to significant discrepancies in the atmospheric 
-!--       surface forcing. For this reason, the child domain
-!--       is initialized with mean soil profiles from the root domain, even if
-!--       no initialization with inifor is set. 
-          init_tsoil_from_parent = .FALSE. 
-          init_msoil_from_parent = .FALSE. 
-          IF ( nested_run )  THEN
-#if defined( __parallel )
-             CALL MPI_ALLREDUCE( init_3d%from_file_tsoil,                      &
-                                 init_tsoil_from_parent,                       &
-                                 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr )
-             CALL MPI_ALLREDUCE( init_3d%from_file_msoil,                      &
-                                 init_msoil_from_parent,                       &
-                                 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr )
-#endif
-          ENDIF
-!
 !--       First, initialize soil temperature and moisture. 
 !--       According to the initialization for surface and soil parameters, 
@@ -4070,110 +4083,248 @@
           ENDDO
 !
-!--       Initialization of soil moisture and temperature from file. 
-!--       In case of no dynamic input file is available for the child domain, 
-!--       transfer soil mean profiles from the root-parent domain onto all 
-!--       child domains.
-          IF ( init_msoil_from_parent )  THEN
-!
-!--          Child domains will be only initialized with horizontally 
-!--          averaged soil profiles in parent domain (for sake of simplicity). 
-!--          If required, average soil data on root parent domain before 
-!--          distribute onto child domains.
-             IF ( init_3d%from_file_msoil  .AND.  init_3d%lod_msoil == 2 )     &
-             THEN
-                ALLOCATE( pr_soil_init(0:init_3d%nzs-1) )
-
-                DO  k = 0, init_3d%nzs-1
-                   pr_soil_init(k) = SUM( init_3d%msoil_3d(k,nys:nyn,nxl:nxr)  )
-                ENDDO
-!
-!--             Allocate 1D array for soil-moisture profile (will not be 
-!--             allocated in lod==2 case). 
-                ALLOCATE( init_3d%msoil_1d(0:init_3d%nzs-1) )
-                init_3d%msoil_1d = 0.0_wp
-#if defined( __parallel )
-                CALL MPI_ALLREDUCE( pr_soil_init(0), init_3d%msoil_1d(0),      &
-                                    SIZE(pr_soil_init),                        &
-                                    MPI_REAL, MPI_SUM, comm2d, ierr )
-#endif
-                init_3d%msoil_1d = init_3d%msoil_1d /                          &
-                                        REAL( ( nx + 1 ) * ( ny + 1), KIND=wp )
-                DEALLOCATE( pr_soil_init )
+!--       Level 2 initialization of the soil, read soil properties from 
+!--       dynamic input file. 
+          IF ( input_pids_dynamic )  THEN
+!
+!--          CPU measurement
+             CALL cpu_log( log_point_s(85), 'NetCDF input init', 'start' )
+#if defined ( __netcdf )
+!
+!--          Open file in read-only mode
+             CALL open_read_file( TRIM( input_file_dynamic ) //                &
+                                  TRIM( coupling_char ), pids_id )
+!
+!--          Inquire all variable names
+             CALL inquire_num_variables( pids_id, num_var_pids )
+!
+!--          Allocate memory to store variable names.
+             ALLOCATE( vars_pids(1:num_var_pids) )
+             CALL inquire_variable_names( pids_id, vars_pids )
+!           
+!--          Read vertical dimension for soil depth. 
+             IF ( check_existence( vars_pids, 'zsoil' ) )                      &
+                CALL get_dimension_length( pids_id, init_3d%nzs, 'zsoil' )
+!           
+!--          Read also the horizontal dimensions required for soil initialization.
+!--          Please note, in case of non-nested runs or in case of root domain, 
+!--          these data is already available, but will be read again for the sake
+!--          of clearness. 
+             CALL get_dimension_length( pids_id, init_3d%nx, 'x'  )
+             CALL get_dimension_length( pids_id, init_3d%ny, 'y'  )
+!           
+!--          Check for correct horizontal and vertical dimension. Please note,
+!--          in case of non-nested runs or in case of root domain, these checks 
+!--          are already performed
+             IF ( init_3d%nx-1 /= nx  .OR.  init_3d%ny-1 /= ny )  THEN
+                message_string = 'Number of horizontal grid points in '//      &
+                                 'dynamic input file does not match ' //       &
+                                 'the number of numeric grid points.'
+                CALL message( 'lsm_init', 'PA0543', 1, 2, 0, 6, 0 )
              ENDIF
-          ENDIF
-          IF ( init_tsoil_from_parent )  THEN
-             IF ( init_3d%from_file_tsoil  .AND.  init_3d%lod_tsoil == 2 )  THEN
-                ALLOCATE( pr_soil_init(0:init_3d%nzs-1) )
-
-                DO  k = 0, init_3d%nzs-1
-                   pr_soil_init(k) = SUM( init_3d%tsoil_3d(k,nys:nyn,nxl:nxr)  )
-                ENDDO
-!
-!--             Allocate 1D array for soil-temperature profile (will not be 
-!--             allocated in lod==2 case). 
-                ALLOCATE( init_3d%tsoil_1d(0:init_3d%nzs-1) )
-                init_3d%tsoil_1d = 0.0_wp
-#if defined( __parallel )
-                CALL MPI_ALLREDUCE( pr_soil_init(0), init_3d%tsoil_1d(0),      &
-                                    SIZE(pr_soil_init),                        &
-                                    MPI_REAL, MPI_SUM, comm2d, ierr )
-#endif
-                init_3d%tsoil_1d = init_3d%tsoil_1d /                          &
-                                        REAL( ( nx + 1 ) * ( ny + 1), KIND=wp )
-                DEALLOCATE( pr_soil_init )
-
+!           
+!--          Read vertical dimensions. Later, these are required for eventual
+!--          inter- and extrapolations of the initialization data.
+             IF ( check_existence( vars_pids, 'zsoil' ) )  THEN
+                ALLOCATE( init_3d%z_soil(1:init_3d%nzs) )
+                CALL get_variable( pids_id, 'zsoil', init_3d%z_soil )
              ENDIF
-          ENDIF
-          IF ( init_msoil_from_parent  .OR.  init_tsoil_from_parent )  THEN
-!
-!--          Distribute soil grid information on file from root to all childs.
-!--          Only process with rank 0 sends the information. 
-#if defined( __parallel )
-             CALL MPI_BCAST( init_3d%nzs,    1,                                &
-                             MPI_INTEGER, 0, MPI_COMM_WORLD, ierr )
-#endif
-
-             IF ( .NOT.  ALLOCATED( init_3d%z_soil ) )                         &
-                ALLOCATE( init_3d%z_soil(1:init_3d%nzs) )
-#if defined( __parallel )
-             CALL MPI_BCAST( init_3d%z_soil, SIZE(init_3d%z_soil),             &
-                             MPI_REAL, 0, MPI_COMM_WORLD, ierr )
-#endif
-          ENDIF
-!
-!--       ALLOCATE arrays on child domains and set control attributes.
-!--       Note, 1d soil profiles are allocated even though soil information
-!--       is already read from dynamic file in one child domain. 
-!--       This case, however, data is not used for further initialization 
-!--       since the LoD=2. 
-          IF ( init_msoil_from_parent )  THEN
-             IF ( .NOT. ALLOCATED( init_3d%msoil_1d ) )  THEN
-                ALLOCATE( init_3d%msoil_1d(0:init_3d%nzs-1) )
-                IF( .NOT. init_3d%from_file_msoil )  init_3d%lod_msoil = 1
+!           
+!--          Read initial data for soil moisture
+             IF ( check_existence( vars_pids, 'init_soil_m' ) )  THEN
+!           
+!--             Read attributes for the fill value and level-of-detail
+                CALL get_attribute( pids_id, char_fill,                        &
+                                    init_3d%fill_msoil,                        &
+                                    .FALSE., 'init_soil_m' )
+                CALL get_attribute( pids_id, char_lod,                         &
+                                    init_3d%lod_msoil,                         &
+                                    .FALSE., 'init_soil_m' )
+!           
+!--             level-of-detail 1 - read initialization profile
+                IF ( init_3d%lod_msoil == 1 )  THEN
+                   ALLOCATE( init_3d%msoil_1d(0:init_3d%nzs-1) )
+            
+                   CALL get_variable( pids_id, 'init_soil_m',                  &
+                                      init_3d%msoil_1d(0:init_3d%nzs-1) )
+!           
+!--             level-of-detail 2 - read 3D initialization data
+                ELSEIF ( init_3d%lod_msoil == 2 )  THEN
+                   ALLOCATE ( init_3d%msoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
+            
+                  CALL get_variable( pids_id, 'init_soil_m',                   &   
+                             init_3d%msoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr),&
+                             nxl, nxr, nys, nyn, 0, init_3d%nzs-1 )
+            
+                ENDIF
                 init_3d%from_file_msoil = .TRUE.
              ENDIF
-          ENDIF
-          IF ( init_tsoil_from_parent )  THEN
-             IF ( .NOT. ALLOCATED( init_3d%tsoil_1d ) )  THEN
-                ALLOCATE( init_3d%tsoil_1d(0:init_3d%nzs-1) )
-                IF( .NOT. init_3d%from_file_tsoil )  init_3d%lod_tsoil = 1
+!           
+!--          Read soil temperature
+             IF ( check_existence( vars_pids, 'init_soil_t' ) )  THEN
+!           
+!--             Read attributes for the fill value and level-of-detail
+                CALL get_attribute( pids_id, char_fill,                        &
+                                    init_3d%fill_tsoil,                        &
+                                    .FALSE., 'init_soil_t' )
+                CALL get_attribute( pids_id, char_lod,                         &
+                                    init_3d%lod_tsoil,                         &
+                                    .FALSE., 'init_soil_t' )
+!           
+!--             level-of-detail 1 - read initialization profile
+                IF ( init_3d%lod_tsoil == 1 )  THEN
+                   ALLOCATE( init_3d%tsoil_1d(0:init_3d%nzs-1) )
+            
+                   CALL get_variable( pids_id, 'init_soil_t',                  &
+                                      init_3d%tsoil_1d(0:init_3d%nzs-1) )
+            
+!           
+!--             level-of-detail 2 - read 3D initialization data
+                ELSEIF ( init_3d%lod_tsoil == 2 )  THEN
+                   ALLOCATE ( init_3d%tsoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
+                   
+                   CALL get_variable( pids_id, 'init_soil_t',                  &   
+                             init_3d%tsoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr),&
+                             nxl, nxr, nys, nyn, 0, init_3d%nzs-1 )
+                ENDIF
                 init_3d%from_file_tsoil = .TRUE.
              ENDIF
-          ENDIF
-!
-!--       Distribute soil profiles from root to all childs
-          IF ( init_msoil_from_parent )  THEN
+!           
+!--          Close input file
+             CALL close_input_file( pids_id )
+#endif      
+!           
+!--          End of CPU measurement
+             CALL cpu_log( log_point_s(85), 'NetCDF input init', 'stop' )
+          ENDIF
+!
+!--       In case no dynamic input is available for a child domain but the 
+!--       parent is initialized with dynamic input file, the different soil 
+!--       states can lead to significant discrepancies in the atmospheric 
+!--       surface forcing. For this reason, the child domain is initialized with
+!--       domain-averaged soil profiles from the root domain, even if
+!--       no initialization with inifor is set. Note, as long as a dynamic 
+!--       input file with soil information is available for the child domain,
+!--       the input file information will be used.
+          IF ( nested_run )  THEN
 #if defined( __parallel )
-             CALL MPI_BCAST( init_3d%msoil_1d, SIZE(init_3d%msoil_1d),         &
+!
+!--          In case of a nested run, first average the soil profiles in the 
+!--          root domain.
+             IF ( pmc_is_rootmodel() )  THEN
+!            
+!--             Child domains will be only initialized with horizontally 
+!--             averaged soil profiles in parent domain (for sake of simplicity). 
+!--             If required, average soil data on root parent domain before the 
+!--             soil profiles are distributed onto the child domains.
+!--             Start with soil moisture.
+                IF ( init_3d%from_file_msoil  .AND.  init_3d%lod_msoil == 2 )  &
+                THEN
+                   ALLOCATE( pr_soil_init(0:init_3d%nzs-1) )
+                   DO  k = 0, init_3d%nzs-1
+                      pr_soil_init(k) = SUM( init_3d%msoil_3d(k,nys:nyn,nxl:nxr)  )
+                   ENDDO
+!            
+!--                Allocate 1D array for soil-moisture profile (will not be 
+!--                allocated in lod==2 case). 
+                   ALLOCATE( init_3d%msoil_1d(0:init_3d%nzs-1) )
+                   init_3d%msoil_1d = 0.0_wp
+                   CALL MPI_ALLREDUCE( pr_soil_init(0), init_3d%msoil_1d(0),   &
+                                       SIZE(pr_soil_init),                     &
+                                       MPI_REAL, MPI_SUM, comm2d, ierr )
+             
+                   init_3d%msoil_1d = init_3d%msoil_1d /                       &
+                                        REAL( ( nx + 1 ) * ( ny + 1), KIND=wp )
+                   DEALLOCATE( pr_soil_init )
+                ENDIF
+! 
+!--             Proceed with soil temperature.
+                IF ( init_3d%from_file_tsoil  .AND.  init_3d%lod_tsoil == 2 )  &
+                THEN
+                   ALLOCATE( pr_soil_init(0:init_3d%nzs-1) )
+             
+                   DO  k = 0, init_3d%nzs-1
+                      pr_soil_init(k) = SUM( init_3d%tsoil_3d(k,nys:nyn,nxl:nxr)  )
+                   ENDDO
+!            
+!--                Allocate 1D array for soil-temperature profile (will not be 
+!--                allocated in lod==2 case). 
+                   ALLOCATE( init_3d%tsoil_1d(0:init_3d%nzs-1) )
+                   init_3d%tsoil_1d = 0.0_wp
+                   CALL MPI_ALLREDUCE( pr_soil_init(0), init_3d%tsoil_1d(0),   &
+                                       SIZE(pr_soil_init),                     &
+                                       MPI_REAL, MPI_SUM, comm2d, ierr )
+                   init_3d%tsoil_1d = init_3d%tsoil_1d /                       &
+                                        REAL( ( nx + 1 ) * ( ny + 1), KIND=wp )
+                   DEALLOCATE( pr_soil_init )
+             
+                ENDIF
+             ENDIF
+!
+!--          Broadcast number of soil layers in root model to all childs. 
+!--          Note, only process 0 in COMM_WORLD is sending. 
+             IF ( pmc_is_rootmodel() )  nzs_root = init_3d%nzs
+             
+             CALL MPI_BCAST( nzs_root, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr )
+!            
+!--          Allocate dummy arrays for soil moisture and temperature profiles 
+!--          on all domains.             
+             ALLOCATE( z_soil_root(1:nzs_root)   )
+             ALLOCATE( m_soil_root(0:nzs_root-1) )
+             ALLOCATE( t_soil_root(0:nzs_root-1) )
+!
+!--          Distribute the mean soil profiles to all child domains.
+             IF ( pmc_is_rootmodel() )  THEN
+                z_soil_root = init_3d%z_soil
+                m_soil_root = init_3d%msoil_1d
+                t_soil_root = init_3d%tsoil_1d
+             ENDIF
+             
+             CALL MPI_BCAST( z_soil_root, SIZE( z_soil_root ),                 &
+                             MPI_REAL, 0, MPI_COMM_WORLD, ierr )               
+             CALL MPI_BCAST( m_soil_root, SIZE( m_soil_root ),                 &
+                             MPI_REAL, 0, MPI_COMM_WORLD, ierr )               
+             CALL MPI_BCAST( t_soil_root, SIZE( t_soil_root ),                 &
                              MPI_REAL, 0, MPI_COMM_WORLD, ierr )
+!
+!--          In the following, the child domains decide whether they take
+!--          the information from the root domain or not. 
+             IF ( .NOT. pmc_is_rootmodel() )  THEN
+!
+!--             If soil moisture or temperature isn't in dynamic input file for
+!--             the child, take the information provided from the root model. 
+!--             Start with z-dimension
+                IF ( .NOT. init_3d%from_file_msoil  .OR.                       &
+                     .NOT. init_3d%from_file_msoil    )  THEN
+                   init_3d%nzs = nzs_root
+                   ALLOCATE( init_3d%z_soil(1:init_3d%nzs) )
+                   init_3d%z_soil(1:init_3d%nzs) = z_soil_root
+                ENDIF
+!               
+!--             Take soil moisture. Note, control flags from_file... and LoD 
+!--             need to be set. 
+                IF ( .NOT. init_3d%from_file_msoil )  THEN
+                   ALLOCATE( init_3d%msoil_1d(0:init_3d%nzs-1) )
+                   init_3d%lod_msoil = 1
+                   init_3d%from_file_msoil = .TRUE.
+                   
+                   init_3d%msoil_1d = m_soil_root             
+                ENDIF
+!               
+!--             Take soil temperature. Note, control flags from_file... and LoD 
+!--             need to be set. 
+                IF (  .NOT. init_3d%from_file_tsoil )  THEN
+                   ALLOCATE( init_3d%tsoil_1d(0:init_3d%nzs-1) )
+                   init_3d%lod_tsoil = 1
+                   init_3d%from_file_tsoil = .TRUE.
+                   
+                   init_3d%tsoil_1d = t_soil_root  
+                ENDIF
+             ENDIF
+             
+             DEALLOCATE( z_soil_root )
+             DEALLOCATE( m_soil_root )
+             DEALLOCATE( t_soil_root )
+          ENDIF
 #endif
-                    
-          ENDIF
-          IF ( init_tsoil_from_parent )  THEN
-#if defined( __parallel )
-             CALL MPI_BCAST( init_3d%tsoil_1d, SIZE(init_3d%tsoil_1d),         &
-                             MPI_REAL, 0, MPI_COMM_WORLD, ierr )
-#endif
-          ENDIF
 !
 !--       Proceed with Level 2 initialization. 
@@ -4185,5 +4336,5 @@
                         surf_lsm_h%pavement_surface(m) )  THEN
 
-                      CALL netcdf_data_input_interpolate(                      &
+                      CALL interpolate_soil_profile(                           &
                                        m_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
                                        init_3d%msoil_1d(:),                    &
@@ -4197,5 +4348,5 @@
                       IF ( surf_lsm_v(l)%vegetation_surface(m)  .OR.           &
                            surf_lsm_v(l)%pavement_surface(m) )  THEN
-                         CALL netcdf_data_input_interpolate(                   &
+                         CALL interpolate_soil_profile(                        &
                                        m_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
                                        init_3d%msoil_1d(:),                    &
@@ -4215,5 +4366,5 @@
 
                       IF ( init_3d%msoil_3d(0,j,i) /= init_3d%fill_msoil )     &
-                         CALL netcdf_data_input_interpolate(                   &
+                         CALL interpolate_soil_profile(                        &
                                        m_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
                                        init_3d%msoil_3d(:,j,i),                &
@@ -4236,5 +4387,5 @@
                          
                          IF ( init_3d%msoil_3d(0,j,i) /= init_3d%fill_msoil )  &
-                            CALL netcdf_data_input_interpolate(                &
+                            CALL interpolate_soil_profile(                     &
                                        m_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
                                        init_3d%msoil_3d(:,j,i),                &
@@ -4255,5 +4406,5 @@
                    IF ( surf_lsm_h%vegetation_surface(m)  .OR.                 &
                         surf_lsm_h%pavement_surface(m) )  THEN
-                      CALL netcdf_data_input_interpolate(                      &
+                      CALL interpolate_soil_profile(                           &
                                        t_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
                                        init_3d%tsoil_1d(:),                    &
@@ -4270,5 +4421,5 @@
                       IF ( surf_lsm_v(l)%vegetation_surface(m)  .OR.           &
                            surf_lsm_v(l)%pavement_surface(m) )  THEN
-                        CALL netcdf_data_input_interpolate(                    &
+                        CALL interpolate_soil_profile(                         &
                                        t_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
                                        init_3d%tsoil_1d(:),                    &
@@ -4292,5 +4443,5 @@
                       
                       IF ( init_3d%tsoil_3d(0,j,i) /= init_3d%fill_tsoil )     &
-                         CALL netcdf_data_input_interpolate(                   &
+                         CALL interpolate_soil_profile(                        &
                                        t_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
                                        init_3d%tsoil_3d(:,j,i),                &
@@ -4316,5 +4467,5 @@
                          
                          IF ( init_3d%tsoil_3d(0,j,i) /= init_3d%fill_tsoil )  &
-                            CALL netcdf_data_input_interpolate(                &
+                            CALL interpolate_soil_profile(                     &
                                        t_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
                                        init_3d%tsoil_3d(:,j,i),                &
@@ -5298,5 +5449,5 @@
 !--                (mass conservation is violated here)
                    DO  k = nzb_soil, nzt_soil
-                      surf_m_soil_p%var_2d(k,m) = MIN( surf_m_soil_p%var_2d(k,m), surf_m_soil_p%var_2d(k,m) )
+                      surf_m_soil_p%var_2d(k,m) = MIN( surf_m_soil_p%var_2d(k,m), surf%m_sat(k,m) )
                       surf_m_soil_p%var_2d(k,m) = MAX( surf_m_soil_p%var_2d(k,m), 0.0_wp )                      
                    ENDDO
@@ -7061,4 +7212,59 @@
 
     END SUBROUTINE calc_z0_water_surface
+    
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!>  Vertical interpolation and extrapolation of 1D soil profile input from 
+!>  dynamic input file onto the numeric vertical soil grid. 
+!------------------------------------------------------------------------------!
+    SUBROUTINE interpolate_soil_profile( var, var_file, z_grid, z_file,        &
+                                         nzb_var, nzt_var, nzb_file, nzt_file )
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  k        !< running index z-direction file
+       INTEGER(iwp) ::  kk       !< running index z-direction stretched model grid
+       INTEGER(iwp) ::  ku       !< upper index bound along z-direction for varialbe from file
+       INTEGER(iwp) ::  nzb_var  !< lower bound of final array
+       INTEGER(iwp) ::  nzt_var  !< upper bound of final array
+       INTEGER(iwp) ::  nzb_file !< lower bound of file array
+       INTEGER(iwp) ::  nzt_file !< upper bound of file array
+
+       REAL(wp), DIMENSION(nzb_var:nzt_var)   ::  z_grid   !< grid levels on numeric grid
+       REAL(wp), DIMENSION(nzb_file:nzt_file) ::  z_file   !< grid levels on file grid
+       REAL(wp), DIMENSION(nzb_var:nzt_var)   ::  var      !< treated variable
+       REAL(wp), DIMENSION(nzb_file:nzt_file) ::  var_file !< temporary variable
+
+       ku = nzt_file
+
+       DO  k = nzb_var, nzt_var
+!
+!--       Determine index on Inifor grid which is closest to the actual height
+          kk = MINLOC( ABS( z_file - z_grid(k) ), DIM = 1 )
+!
+!--       If closest index on Inifor grid is smaller than top index,
+!--       interpolate the data
+          IF ( kk < nzt_file )  THEN
+             IF ( z_file(kk) - z_grid(k) <= 0.0_wp )  THEN
+                var(k) = var_file(kk) + ( var_file(kk+1) - var_file(kk) ) /    &
+                                        ( z_file(kk+1)   - z_file(kk)   ) *    &
+                                        ( z_grid(k)      - z_file(kk)   )
+
+             ELSEIF ( z_file(kk) - z_grid(k) > 0.0_wp )  THEN
+                var(k) = var_file(kk-1) + ( var_file(kk) - var_file(kk-1) ) /  &
+                                          ( z_file(kk)   - z_file(kk-1)   ) *  &
+                                          ( z_grid(k)    - z_file(kk-1)   )
+             ENDIF
+!
+!--       Extrapolate if actual height is above the highest Inifor level by the last value
+          ELSE
+             var(k) = var_file(ku)
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE interpolate_soil_profile
 
 !
Index: palm/trunk/SOURCE/netcdf_data_input_mod.f90
===================================================================
--- palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 4256)
+++ palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 4258)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! - Migrate input of soil temperature and moisture to land-surface model.
+! - Remove interpolate routines and move the only required subroutine to 
+!   land-surface model.
+! 
+! 4247 2019-09-30 10:18:24Z pavelkrc
 ! Add reading and processing of building_surface_pars
 ! 
@@ -608,11 +613,4 @@
     PRIVATE
 
-    INTERFACE netcdf_data_input_interpolate
-       MODULE PROCEDURE netcdf_data_input_interpolate_1d
-       MODULE PROCEDURE netcdf_data_input_interpolate_1d_soil
-       MODULE PROCEDURE netcdf_data_input_interpolate_2d
-       MODULE PROCEDURE netcdf_data_input_interpolate_3d
-    END INTERFACE netcdf_data_input_interpolate
-
     INTERFACE netcdf_data_input_check_dynamic
        MODULE PROCEDURE netcdf_data_input_check_dynamic
@@ -650,8 +648,4 @@
     END INTERFACE netcdf_data_input_init_3d
     
-    INTERFACE netcdf_data_input_init_lsm
-       MODULE PROCEDURE netcdf_data_input_init_lsm
-    END INTERFACE netcdf_data_input_init_lsm
-
     INTERFACE netcdf_data_input_surface_data
        MODULE PROCEDURE netcdf_data_input_surface_data
@@ -733,12 +727,14 @@
 !
 !-- Public subroutines
-    PUBLIC netcdf_data_input_check_dynamic, netcdf_data_input_check_static,    &
+    PUBLIC netcdf_data_input_check_dynamic,                                    &
+           netcdf_data_input_check_static,                                     &
            netcdf_data_input_chemistry_data,                                   &
            get_dimension_length,                                               &
            netcdf_data_input_inquire_file,                                     &
-           netcdf_data_input_init, netcdf_data_input_init_lsm,                 &
-           netcdf_data_input_init_3d, netcdf_data_input_att,                   &
-           netcdf_data_input_interpolate,                                      &
-           netcdf_data_input_surface_data, netcdf_data_input_topo,             &
+           netcdf_data_input_init,                                             &
+           netcdf_data_input_init_3d,                                          &
+           netcdf_data_input_att,                                              &
+           netcdf_data_input_surface_data,                                     &
+           netcdf_data_input_topo,                                             &
            netcdf_data_input_var,                                              &
            get_attribute,                                                      &
@@ -2914,14 +2910,14 @@
        IF ( init_3d%nx-1 /= nx  .OR.  init_3d%nxu-1 /= nx - 1  .OR.            &
             init_3d%ny-1 /= ny  .OR.  init_3d%nyv-1 /= ny - 1 )  THEN
-          message_string = 'Number of inifor horizontal grid points  '//       &
-                           'does not match the number of numeric grid '//      &
-                           'points.'
+          message_string = 'Number of horizontal grid points in '//            &
+                           'dynamic input file does not match ' //             &
+                           'the number of numeric grid points.'
           CALL message( 'netcdf_data_input_mod', 'PA0543', 1, 2, 0, 6, 0 )
        ENDIF
 
        IF ( init_3d%nzu /= nz )  THEN
-          message_string = 'Number of inifor vertical grid points ' //         &
-                           'does not match the number of numeric grid '//      &
-                           'points.'
+          message_string = 'Number of vertical grid points in '//              &
+                           'dynamic input file does not match ' //             &
+                           'the number of numeric grid points.'
           CALL message( 'netcdf_data_input_mod', 'PA0543', 1, 2, 0, 6, 0 )
        ENDIF
@@ -3353,143 +3349,4 @@
 
     END SUBROUTINE netcdf_data_input_init_3d
-    
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Reads initialization data of u, v, w, pt, q, geostrophic wind components,
-!> as well as soil moisture and soil temperature, derived from larger-scale
-!> model (COSMO) by Inifor.
-!------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_init_lsm
-
-       USE control_parameters,                                                 &
-           ONLY:  message_string
-
-       USE indices,                                                            &
-           ONLY:  nx, nxl, nxr, ny, nyn, nys
-
-       IMPLICIT NONE
-
-       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names !< string containing all variables on file
-      
-       INTEGER(iwp) ::  id_dynamic !< NetCDF id of dynamic input file
-       INTEGER(iwp) ::  num_vars   !< number of variables in netcdf input file
-
-!
-!--    Skip routine if no input file with dynamic input data is available.
-       IF ( .NOT. input_pids_dynamic )  RETURN
-!
-!--    CPU measurement
-       CALL cpu_log( log_point_s(85), 'NetCDF input init', 'start' )
-
-#if defined ( __netcdf )
-!
-!--    Open file in read-only mode
-       CALL open_read_file( TRIM( input_file_dynamic ) //                      &
-                            TRIM( coupling_char ), id_dynamic )
-
-!
-!--    At first, inquire all variable names.
-       CALL inquire_num_variables( id_dynamic, num_vars )
-!
-!--    Allocate memory to store variable names.
-       ALLOCATE( var_names(1:num_vars) )
-       CALL inquire_variable_names( id_dynamic, var_names )
-!
-!--    Read vertical dimension for soil depth. 
-       IF ( check_existence( var_names, 'zsoil' ) )                            &
-          CALL get_dimension_length( id_dynamic, init_3d%nzs, 'zsoil' )
-!
-!--    Read also the horizontal dimensions required for soil initialization.
-!--    Please note, in case of non-nested runs or in case of root domain, 
-!--    these data is already available, but will be read again for the sake
-!--    of clearness. 
-       CALL get_dimension_length( id_dynamic, init_3d%nx, 'x'  )
-       CALL get_dimension_length( id_dynamic, init_3d%ny, 'y'  )
-!
-!--    Check for correct horizontal and vertical dimension. Please note,
-!--    in case of non-nested runs or in case of root domain, these checks 
-!--    are already performed
-       IF ( init_3d%nx-1 /= nx  .OR.  init_3d%ny-1 /= ny )  THEN
-          message_string = 'Number of inifor horizontal grid points  '//       &
-                           'does not match the number of numeric grid points.'
-          CALL message( 'netcdf_data_input_mod', 'PA0543', 1, 2, 0, 6, 0 )
-       ENDIF
-!
-!--    Read vertical dimensions. Later, these are required for eventual
-!--    inter- and extrapolations of the initialization data.
-       IF ( check_existence( var_names, 'zsoil' ) )  THEN
-          ALLOCATE( init_3d%z_soil(1:init_3d%nzs) )
-          CALL get_variable( id_dynamic, 'zsoil', init_3d%z_soil )
-       ENDIF
-!
-!--    Read initial data for soil moisture
-       IF ( check_existence( var_names, 'init_soil_m' ) )  THEN
-!
-!--       Read attributes for the fill value and level-of-detail
-          CALL get_attribute( id_dynamic, char_fill,                           &
-                              init_3d%fill_msoil,                              &
-                              .FALSE., 'init_soil_m' )
-          CALL get_attribute( id_dynamic, char_lod,                            &
-                              init_3d%lod_msoil,                               &
-                              .FALSE., 'init_soil_m' )
-!
-!--       level-of-detail 1 - read initialization profile
-          IF ( init_3d%lod_msoil == 1 )  THEN
-             ALLOCATE( init_3d%msoil_1d(0:init_3d%nzs-1) )
-
-             CALL get_variable( id_dynamic, 'init_soil_m',                     &
-                                init_3d%msoil_1d(0:init_3d%nzs-1) )
-!
-!--       level-of-detail 2 - read 3D initialization data
-          ELSEIF ( init_3d%lod_msoil == 2 )  THEN
-             ALLOCATE ( init_3d%msoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
-
-            CALL get_variable( id_dynamic, 'init_soil_m',                      &   
-                             init_3d%msoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr),&
-                             nxl, nxr, nys, nyn, 0, init_3d%nzs-1 )
-
-          ENDIF
-          init_3d%from_file_msoil = .TRUE.
-       ENDIF
-!
-!--    Read soil temperature
-       IF ( check_existence( var_names, 'init_soil_t' ) )  THEN
-!
-!--       Read attributes for the fill value and level-of-detail
-          CALL get_attribute( id_dynamic, char_fill,                           &
-                              init_3d%fill_tsoil,                              &
-                              .FALSE., 'init_soil_t' )
-          CALL get_attribute( id_dynamic, char_lod,                            &
-                              init_3d%lod_tsoil,                               &
-                              .FALSE., 'init_soil_t' )
-!
-!--       level-of-detail 1 - read initialization profile
-          IF ( init_3d%lod_tsoil == 1 )  THEN
-             ALLOCATE( init_3d%tsoil_1d(0:init_3d%nzs-1) )
-
-             CALL get_variable( id_dynamic, 'init_soil_t',                     &
-                                init_3d%tsoil_1d(0:init_3d%nzs-1) )
-
-!
-!--       level-of-detail 2 - read 3D initialization data
-          ELSEIF ( init_3d%lod_tsoil == 2 )  THEN
-             ALLOCATE ( init_3d%tsoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
-             
-             CALL get_variable( id_dynamic, 'init_soil_t',                     &   
-                             init_3d%tsoil_3d(0:init_3d%nzs-1,nys:nyn,nxl:nxr),&
-                             nxl, nxr, nys, nyn, 0, init_3d%nzs-1 )
-          ENDIF
-          init_3d%from_file_tsoil = .TRUE.
-       ENDIF
-!
-!--    Close input file
-       CALL close_input_file( id_dynamic )
-#endif
-!
-!--    End of CPU measurement
-       CALL cpu_log( log_point_s(85), 'NetCDF input init', 'stop' )
-
-    END SUBROUTINE netcdf_data_input_init_lsm    
 
 !------------------------------------------------------------------------------!
@@ -3584,5 +3441,5 @@
                 message_string = 'Reading 3D building data - vertical ' //     &
                                  'coordinate do not match numeric grid.'
-                CALL message( 'netcdf_data_input_mod', 'PA0553', 2, 2, myid, 6, 0 )
+                CALL message( 'netcdf_data_input_mod', 'PA0553', 2, 2, 0, 6, 0 )
              ENDIF
           ENDIF
@@ -4227,262 +4084,4 @@
 
     END SUBROUTINE resize_array_4d_real
-    
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Vertical interpolation and extrapolation of 1D variables.
-!------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_interpolate_1d( var, z_grid, z_file)
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  k       !< running index z-direction file
-       INTEGER(iwp) ::  kk      !< running index z-direction stretched model grid
-       INTEGER(iwp) ::  kl      !< lower index bound along z-direction
-       INTEGER(iwp) ::  ku      !< upper index bound along z-direction
-
-       REAL(wp), DIMENSION(:) ::  z_grid                  !< grid levels on numeric grid
-       REAL(wp), DIMENSION(:) ::  z_file                  !< grid levels on file grid
-       REAL(wp), DIMENSION(:), INTENT(INOUT) ::  var      !< treated variable
-       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var_tmp  !< temporary variable
-
-
-       kl = LBOUND(var,1)
-       ku = UBOUND(var,1)
-       ALLOCATE( var_tmp(kl:ku) )
-
-       DO  k = kl, ku
-
-          kk = MINLOC( ABS( z_file - z_grid(k) ), DIM = 1 )
-
-          IF ( kk < ku )  THEN
-             IF ( z_file(kk) - z_grid(k) <= 0.0_wp )  THEN
-                var_tmp(k) = var(kk) +                                         &
-                                       ( var(kk+1)        - var(kk)    ) /     &
-                                       ( z_file(kk+1)     - z_file(kk) ) *     &
-                                       ( z_grid(k)        - z_file(kk) )
-
-             ELSEIF ( z_file(kk) - z_grid(k) > 0.0_wp )  THEN
-                var_tmp(k) = var(kk-1) +                                       &
-                                         ( var(kk)     - var(kk-1)    ) /      &
-                                         ( z_file(kk)  - z_file(kk-1) ) *      &
-                                         ( z_grid(k)   - z_file(kk-1) )
-             ENDIF
-!
-!--       Extrapolate
-          ELSE
-
-             var_tmp(k) = var(ku) +   ( var(ku)    - var(ku-1)      ) /        &
-                                      ( z_file(ku) - z_file(ku-1)   ) *        &
-                                      ( z_grid(k)  - z_file(ku)     )
-
-          ENDIF
-
-       ENDDO
-       var(:) = var_tmp(:)
-
-       DEALLOCATE( var_tmp )
-
-
-    END SUBROUTINE netcdf_data_input_interpolate_1d
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Vertical interpolation and extrapolation of 1D variables from Inifor grid
-!> onto Palm grid, where both have same dimension. Please note, the passed
-!> paramter list in 1D version is different compared to 2D version.
-!------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_interpolate_1d_soil( var, var_file,           &
-                                                      z_grid, z_file,          &
-                                                      nzb_var, nzt_var,        &
-                                                      nzb_file, nzt_file )
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  k        !< running index z-direction file
-       INTEGER(iwp) ::  kk       !< running index z-direction stretched model grid
-       INTEGER(iwp) ::  ku       !< upper index bound along z-direction for varialbe from file
-       INTEGER(iwp) ::  nzb_var  !< lower bound of final array
-       INTEGER(iwp) ::  nzt_var  !< upper bound of final array
-       INTEGER(iwp) ::  nzb_file !< lower bound of file array
-       INTEGER(iwp) ::  nzt_file !< upper bound of file array
-
-!        LOGICAL, OPTIONAL ::  zsoil !< flag indicating reverse z-axis, i.e. zsoil instead of height, e.g. in case of ocean or soil
-
-       REAL(wp), DIMENSION(nzb_var:nzt_var)   ::  z_grid   !< grid levels on numeric grid
-       REAL(wp), DIMENSION(nzb_file:nzt_file) ::  z_file   !< grid levels on file grid
-       REAL(wp), DIMENSION(nzb_var:nzt_var)   ::  var      !< treated variable
-       REAL(wp), DIMENSION(nzb_file:nzt_file) ::  var_file !< temporary variable
-
-       ku = nzt_file
-
-       DO  k = nzb_var, nzt_var
-!
-!--       Determine index on Inifor grid which is closest to the actual height
-          kk = MINLOC( ABS( z_file - z_grid(k) ), DIM = 1 )
-!
-!--       If closest index on Inifor grid is smaller than top index,
-!--       interpolate the data
-          IF ( kk < nzt_file )  THEN
-             IF ( z_file(kk) - z_grid(k) <= 0.0_wp )  THEN
-                var(k) = var_file(kk) + ( var_file(kk+1) - var_file(kk) ) /    &
-                                        ( z_file(kk+1)   - z_file(kk)   ) *    &
-                                        ( z_grid(k)      - z_file(kk)   )
-
-             ELSEIF ( z_file(kk) - z_grid(k) > 0.0_wp )  THEN
-                var(k) = var_file(kk-1) + ( var_file(kk) - var_file(kk-1) ) /  &
-                                          ( z_file(kk)   - z_file(kk-1)   ) *  &
-                                          ( z_grid(k)    - z_file(kk-1)   )
-             ENDIF
-!
-!--       Extrapolate if actual height is above the highest Inifor level by the last value
-          ELSE
-             var(k) = var_file(ku)
-          ENDIF
-
-       ENDDO
-
-    END SUBROUTINE netcdf_data_input_interpolate_1d_soil
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Vertical interpolation and extrapolation of 2D variables at lateral boundaries.
-!------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_interpolate_2d( var, z_grid, z_file)
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i       !< running index x- or y -direction
-       INTEGER(iwp) ::  il      !< lower index bound along x- or y-direction
-       INTEGER(iwp) ::  iu      !< upper index bound along x- or y-direction
-       INTEGER(iwp) ::  k       !< running index z-direction file
-       INTEGER(iwp) ::  kk      !< running index z-direction stretched model grid
-       INTEGER(iwp) ::  kl      !< lower index bound along z-direction
-       INTEGER(iwp) ::  ku      !< upper index bound along z-direction
-
-       REAL(wp), DIMENSION(:) ::  z_grid                  !< grid levels on numeric grid
-       REAL(wp), DIMENSION(:) ::  z_file                  !< grid levels on file grid
-       REAL(wp), DIMENSION(:,:), INTENT(INOUT) ::  var    !< treated variable
-       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var_tmp  !< temporary variable
-
-
-       il = LBOUND(var,2)
-       iu = UBOUND(var,2)
-       kl = LBOUND(var,1)
-       ku = UBOUND(var,1)
-       ALLOCATE( var_tmp(kl:ku) )
-
-       DO  i = il, iu
-          DO  k = kl, ku
-
-             kk = MINLOC( ABS( z_file - z_grid(k) ), DIM = 1 )
-
-             IF ( kk < ku )  THEN
-                IF ( z_file(kk) - z_grid(k) <= 0.0_wp )  THEN
-                   var_tmp(k) = var(kk,i) +                                    &
-                                          ( var(kk+1,i)      - var(kk,i)  ) /  &
-                                          ( z_file(kk+1)     - z_file(kk) ) *  &
-                                          ( z_grid(k)        - z_file(kk) )
-
-                ELSEIF ( z_file(kk) - z_grid(k) > 0.0_wp )  THEN
-                   var_tmp(k) = var(kk-1,i) +                                  &
-                                            ( var(kk,i)   - var(kk-1,i)  ) /   &
-                                            ( z_file(kk)  - z_file(kk-1) ) *   &
-                                            ( z_grid(k)   - z_file(kk-1) )
-                ENDIF
-!
-!--          Extrapolate
-             ELSE
-
-                var_tmp(k) = var(ku,i) + ( var(ku,i)  - var(ku-1,i)    ) /     &
-                                         ( z_file(ku) - z_file(ku-1)   ) *     &
-                                         ( z_grid(k)  - z_file(ku)     )
-
-             ENDIF
-
-          ENDDO
-          var(:,i) = var_tmp(:)
-
-       ENDDO
-
-       DEALLOCATE( var_tmp )
-
-
-    END SUBROUTINE netcdf_data_input_interpolate_2d
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Vertical interpolation and extrapolation of 3D variables.
-!------------------------------------------------------------------------------!
-    SUBROUTINE netcdf_data_input_interpolate_3d( var, z_grid, z_file )
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i       !< running index x-direction
-       INTEGER(iwp) ::  il      !< lower index bound along x-direction
-       INTEGER(iwp) ::  iu      !< upper index bound along x-direction
-       INTEGER(iwp) ::  j       !< running index y-direction
-       INTEGER(iwp) ::  jl      !< lower index bound along x-direction
-       INTEGER(iwp) ::  ju      !< upper index bound along x-direction
-       INTEGER(iwp) ::  k       !< running index z-direction file
-       INTEGER(iwp) ::  kk      !< running index z-direction stretched model grid
-       INTEGER(iwp) ::  kl      !< lower index bound along z-direction
-       INTEGER(iwp) ::  ku      !< upper index bound along z-direction
-
-       REAL(wp), DIMENSION(:) ::  z_grid                      !< grid levels on numeric grid
-       REAL(wp), DIMENSION(:) ::  z_file                      !< grid levels on file grid
-       REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) ::  var      !< treated variable
-       REAL(wp), DIMENSION(:), ALLOCATABLE       ::  var_tmp  !< temporary variable
-
-       il = LBOUND(var,3)
-       iu = UBOUND(var,3)
-       jl = LBOUND(var,2)
-       ju = UBOUND(var,2)
-       kl = LBOUND(var,1)
-       ku = UBOUND(var,1)
-
-       ALLOCATE( var_tmp(kl:ku) )
-
-       DO  i = il, iu
-          DO  j = jl, ju
-             DO  k = kl, ku
-
-                kk = MINLOC( ABS( z_file - z_grid(k) ), DIM = 1 )
-
-                IF ( kk < ku )  THEN
-                   IF ( z_file(kk) - z_grid(k) <= 0.0_wp )  THEN
-                      var_tmp(k) = var(kk,j,i) +                               &
-                                             ( var(kk+1,j,i) - var(kk,j,i) ) / &
-                                             ( z_file(kk+1)  - z_file(kk)  ) * &
-                                             ( z_grid(k)     - z_file(kk)  )
-
-                   ELSEIF ( z_file(kk) - z_grid(k) > 0.0_wp )  THEN
-                      var_tmp(k) = var(kk-1,j,i) +                             &
-                                             ( var(kk,j,i) - var(kk-1,j,i) ) / &
-                                             ( z_file(kk)  - z_file(kk-1)  ) * &
-                                             ( z_grid(k)   - z_file(kk-1)  )
-                   ENDIF
-!
-!--             Extrapolate
-                ELSE
-                   var_tmp(k) = var(ku,j,i) +                                  &
-                                       ( var(ku,j,i)  - var(ku-1,j,i)   ) /    &
-                                       ( z_file(ku)   - z_file(ku-1)    ) *    &
-                                       ( z_grid(k)    - z_file(ku)      )
-
-                ENDIF
-             ENDDO
-             var(:,j,i) = var_tmp(:)
-          ENDDO
-       ENDDO
-
-       DEALLOCATE( var_tmp )
-
-
-    END SUBROUTINE netcdf_data_input_interpolate_3d
 
 !------------------------------------------------------------------------------!
Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4256)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4258)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! Check if any LAD is prescribed when plant-canopy model is applied.
+! 
+! 4226 2019-09-10 17:03:24Z suehring
 ! Bugfix, missing initialization of heating rate
 ! 
@@ -941,7 +944,8 @@
        INTEGER(iwp) ::  m   !< running index
 
-       REAL(wp) ::  int_bpdf        !< vertical integral for lad-profile construction
+       REAL(wp) ::  canopy_height   !< canopy height for lad-profile construction
        REAL(wp) ::  gradient        !< gradient for lad-profile construction
-       REAL(wp) ::  canopy_height   !< canopy height for lad-profile construction
+       REAL(wp) ::  int_bpdf        !< vertical integral for lad-profile construction      
+       REAL(wp) ::  lad_max         !< maximum LAD value in the model domain, used to perform a check
 
        IF ( debug_output )  CALL debug_message( 'pcm_init', 'start' )
@@ -1115,4 +1119,18 @@
  
        END SELECT
+!
+!--    Check that at least one grid point has an LAD /= 0, else this may 
+!--    cause errors in the radiation model. 
+       lad_max = MAXVAL( lad_s )
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( MPI_IN_PLACE, lad_max, 1, MPI_REAL, MPI_MAX,        &
+                           comm2d, ierr)
+#endif
+       IF ( lad_max <= 0.0_wp )  THEN
+          message_string = 'Plant-canopy model is switched-on but no ' //      &
+                           'plant canopy is present in the model domain.'
+          CALL message( 'pcm_init', 'PA0685', 1, 2, 0, 6, 0 )
+       ENDIF
+    
 !
 !--    Initialize 2D index array indicating canopy top index. 
@@ -1301,5 +1319,5 @@
 
 !
-!--    Try to find radiation model package
+!--    Try to find plant-canopy model package
        REWIND ( 11 )
        line = ' '
@@ -1314,5 +1332,5 @@
 
 !
-!--    Set flag that indicates that the radiation model is switched on
+!--    Set flag that indicates that the plant-canopy model is switched on
        plant_canopy = .TRUE.
 
@@ -1341,5 +1359,5 @@
 
 !
-!--    Set flag that indicates that the radiation model is switched on
+!--    Set flag that indicates that the plant-canopy model is switched on
        plant_canopy = .TRUE.
 
Index: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 4256)
+++ palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 4258)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! Initialization of Obukhov lenght also at vertical surfaces (if allocated).
+! 
+! 4237 2019-09-25 11:33:42Z knoop
 ! Added missing OpenMP directives
 ! 
@@ -554,7 +557,5 @@
 ! Description:
 ! ------------
-!> Initializing actions for the surface layer routine. Basically, this involves
-!> the preparation of a lookup table for the the bulk Richardson number vs 
-!> Obukhov length L when using the lookup table method.
+!> Initializing actions for the surface layer routine.
 !------------------------------------------------------------------------------!
     SUBROUTINE init_surface_layer_fluxes
@@ -562,4 +563,5 @@
        IMPLICIT NONE
 
+       INTEGER(iwp) ::  l  !< running index for vertical surface orientation
 
        CALL location_message( 'initializing surface layer', 'start' )
@@ -572,4 +574,14 @@
           IF ( surf_lsm_h%ns    >= 1 )  surf_lsm_h%ol    = 1.0E10_wp
           IF ( surf_usm_h%ns    >= 1 )  surf_usm_h%ol    = 1.0E10_wp
+          
+          DO  l = 0, 3
+             IF ( surf_def_v(l)%ns >= 1  .AND.                                 &
+                  ALLOCATED( surf_def_v(l)%ol ) )  surf_def_v(l)%ol = 1.0E10_wp
+             IF ( surf_lsm_v(l)%ns >= 1  .AND.                                 &
+                  ALLOCATED( surf_lsm_v(l)%ol ) )  surf_lsm_v(l)%ol = 1.0E10_wp
+             IF ( surf_usm_v(l)%ns >= 1  .AND.                                 &
+                  ALLOCATED( surf_usm_v(l)%ol ) )  surf_usm_v(l)%ol = 1.0E10_wp  
+          ENDDO
+          
        ENDIF
 
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4256)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4258)
@@ -28,4 +28,9 @@
 ! -----------------
 ! $Id$
+! - Add checks to ensure that relative fractions of walls, windowns and green
+!   surfaces sum-up to one.
+! - Revise message calls dealing with local checks.
+! 
+! 4245 2019-09-30 08:40:37Z pavelkrc
 ! Initialize explicit per-surface parameters from building_surface_pars
 ! 
@@ -2355,6 +2360,6 @@
                 THEN
                    WRITE( message_string, * ) 'building_type = is out of ' //  &
-                                        'the valid range at (j,i) = ', j, i
-                   CALL message( 'usm_check_parameters', 'PA0529', 2, 2, 0, 6, 0 )
+                                              'the valid range at (j,i) = ', j, i
+                   CALL message( 'usm_check_parameters', 'PA0529', 2, 2, myid, 6, 0 )
                 ENDIF
              ENDDO
@@ -4434,5 +4439,7 @@
                                     building_pars_f%pars_xy(ind_lambda_surf,j,i)
               
+              write(9,*) m, SUM( surf_usm_h%frac(:,m) ), "indiv", surf_usm_h%frac(0,m), surf_usm_h%frac(1,m), surf_usm_h%frac(2,m) 
            ENDDO
+           flush(9)
 
 
@@ -5028,4 +5035,29 @@
            ENDDO
         ENDIF
+!
+!--     Run further checks to ensure that the respecitve material fractions are 
+!--     prescribed properly.
+        DO  m = 1, surf_usm_h%ns
+           IF ( SUM( surf_usm_h%frac(:,m) ) /= 1.0_wp )  THEN
+              WRITE(message_string,*) 'The relative material fractions do ' // &
+                                      'not sum-up to one at horizotal ' //     &
+                                      'surface. (i,j) = ',                     &
+                                      surf_usm_h%i(m), surf_usm_h%j(m)
+              CALL message( 'urban_surface_model_mod', 'PA0686', 2, 2, myid, 6, 0 )
+           ENDIF
+        ENDDO
+        
+        DO  l = 0, 3
+           DO  m = 1, surf_usm_v(l)%ns
+              IF ( SUM( surf_usm_v(l)%frac(:,m) ) /= 1.0_wp )  THEN
+                 WRITE(message_string,*)                                       &
+                                      'The relative material fractions do ' // &
+                                      'not sum-up to one at vertical ' //      &
+                                      'surface. (i,j) = ',                     &
+                                      surf_usm_v(l)%i(m), surf_usm_v(l)%j(m)
+                 CALL message( 'urban_surface_model_mod', 'PA0686', 2, 2, myid, 6, 0 )
+              ENDIF
+           ENDDO
+        ENDDO
 !       
 !--     Read the surface_types array.
@@ -5059,5 +5091,5 @@
                             surf_usm_h%i(m), surf_usm_h%j(m)
              CALL message( 'urban_surface_model_mod', 'PA0503',                &
-                            0, 0, 0, 6, 0 )
+                            0, 0, myid, 6, 0 )
           ENDIF
           IF ( surf_usm_h%z0h(m) >= surf_usm_h%z_mo(m) )  THEN
@@ -5071,5 +5103,5 @@
                             surf_usm_h%i(m), surf_usm_h%j(m)
              CALL message( 'urban_surface_model_mod', 'PA0507',                &
-                            0, 0, 0, 6, 0 )
+                            0, 0, myid, 6, 0 )
           ENDIF          
        ENDDO
@@ -5087,5 +5119,5 @@
                             surf_usm_v(l)%j(m)+surf_usm_v(l)%joff
                 CALL message( 'urban_surface_model_mod', 'PA0503',              &
-                            0, 0, 0, 6, 0 )
+                            0, 0, myid, 6, 0 )
              ENDIF
              IF ( surf_usm_v(l)%z0h(m) >= surf_usm_v(l)%z_mo(m) )  THEN
@@ -5100,5 +5132,5 @@
                             surf_usm_v(l)%j(m)+surf_usm_v(l)%joff
                 CALL message( 'urban_surface_model_mod', 'PA0507',               &
-                            0, 0, 0, 6, 0 )
+                            0, 0, myid, 6, 0 )
              ENDIF
           ENDDO
