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