Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 3050)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 3051)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Move initialization call for nudging and 1D/3D offline nesting. 
+! Revise initialization with inifor data.
+! 
+! 3045 2018-05-28 07:55:41Z Giersch
 ! Error messages revised
 ! 
@@ -1116,16 +1120,4 @@
 
 !
-!-- Initialize nudging if required
-    IF ( nudging )  THEN
-       CALL nudge_init
-    ENDIF
-
-!
-!-- Initialize reading of large scale forcing from external file - if required
-    IF ( large_scale_forcing  .OR.  forcing )  THEN
-       CALL lsf_init
-    ENDIF
-
-!
 !-- Allocate arrays containing the RK coefficient for calculation of 
 !-- perturbation pressure and turbulent fluxes. At this point values are
@@ -1166,5 +1158,6 @@
           CALL location_message( 'initializing with INIFOR', .FALSE. )
 !
-!--       Read initial 1D profiles from NetCDF file if available. 
+!--       Read initial 1D profiles or 3D data from NetCDF file, depending 
+!--       on the provided level-of-detail. 
 !--       At the moment, only u, v, w, pt and q are provided. 
           CALL netcdf_data_input_init_3d
@@ -1177,19 +1170,27 @@
           IF ( ANY( zu(1:nzt+1) /= init_3d%zu_atmos(1:init_3d%nzu) ) )  THEN
 
-             CALL netcdf_data_input_interpolate( init_3d%u_init(nzb+1:nzt+1),  &
+             IF( init_3d%lod_u == 1 )                                          &
+                CALL netcdf_data_input_interpolate(                            &
+                                                 init_3d%u_init(nzb+1:nzt+1),  &
                                                  zu(nzb+1:nzt+1),              &
                                                  init_3d%zu_atmos )
-             CALL netcdf_data_input_interpolate( init_3d%v_init(nzb+1:nzt+1),  &
+                                                 
+             IF( init_3d%lod_v == 1 )                                          &
+                CALL netcdf_data_input_interpolate(                            &
+                                                 init_3d%v_init(nzb+1:nzt+1),  &
                                                  zu(nzb+1:nzt+1),              &
                                                  init_3d%zu_atmos )
+                                                 
 !              CALL netcdf_data_input_interpolate( init_3d%w_init(nzb+1:nzt),    &
 !                                                  zw(nzb+1:nzt),                &
 !                                                  init_3d%zw_atmos )
-             IF ( .NOT. neutral )                                              &
+
+             IF ( .NOT. neutral  .AND.  init_3d%lod_pt == 1 )                  &
                 CALL netcdf_data_input_interpolate(                            &
                                              init_3d%pt_init(nzb+1:nzt+1),     &
                                              zu(nzb+1:nzt+1),                  &
                                              init_3d%zu_atmos )
-             IF ( humidity )                                                   &
+                                             
+             IF ( humidity  .AND.  init_3d%lod_q == 1 )                        &
                 CALL netcdf_data_input_interpolate(                            &
                                              init_3d%q_init(nzb+1:nzt+1),      &
@@ -1197,9 +1198,12 @@
                                              init_3d%zu_atmos )
           ENDIF
-
-          u_init = init_3d%u_init
-          v_init = init_3d%v_init    
-          IF( .NOT. neutral )  pt_init = init_3d%pt_init
-          IF( humidity      )  q_init  = init_3d%q_init
+!
+!--       In case of LOD=1, initialize 1D profiles and 3D data.  
+          IF( init_3d%lod_u == 1 )  u_init = init_3d%u_init
+          IF( init_3d%lod_v == 1 )  v_init = init_3d%v_init    
+          IF( .NOT. neutral  .AND.  init_3d%lod_pt == 1 )                      &
+             pt_init = init_3d%pt_init
+          IF( humidity  .AND.  init_3d%lod_q == 1 )                            &
+             q_init  = init_3d%q_init
 
 !
@@ -1210,8 +1214,10 @@
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                u(:,j,i) = u_init(:)
-                v(:,j,i) = v_init(:)
-                IF( .NOT. neutral )  pt(:,j,i) = pt_init(:)
-                IF( humidity      )  q(:,j,i)  = q_init(:)
+                IF( init_3d%lod_u == 1 )  u(:,j,i) = u_init(:)
+                IF( init_3d%lod_v == 1 )  v(:,j,i) = v_init(:)
+                IF( .NOT. neutral  .AND.  init_3d%lod_pt == 1 )                &
+                   pt(:,j,i) = pt_init(:)
+                IF( humidity  .AND.  init_3d%lod_q == 1 )                      &
+                   q(:,j,i)  = q_init(:)
              ENDDO
           ENDDO
@@ -2287,4 +2293,13 @@
 
 !
+!-- Initialize nudging if required
+    IF ( nudging )  CALL nudge_init
+
+!
+!-- Initialize 1D/3D offline-nesting with COSMO model and read data from 
+!-- external file.
+    IF ( large_scale_forcing  .OR.  forcing )  CALL lsf_init
+    
+!
 !-- Initialize surface forcing corresponding to large-scale forcing. Therein, 
 !-- initialize heat-fluxes, etc. via datatype. Revise it later!
Index: /palm/trunk/SOURCE/init_grid.f90
===================================================================
--- /palm/trunk/SOURCE/init_grid.f90	(revision 3050)
+++ /palm/trunk/SOURCE/init_grid.f90	(revision 3051)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Minor bugfix concerning mapping 3D buildings on top of terrain
+! 
+! 3045 2018-05-28 07:55:41Z Giersch
 ! Error messages revised
 ! 
@@ -926,4 +929,8 @@
     topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
 !
+!-- In order to map topography on PALM grid also in case of ocean simulations,
+!-- pre-calculate an offset value.
+    ocean_offset = MERGE( zw(0), 0.0_wp, ocean )
+!
 !-- Reference buildings on top of orography. This is not necessary
 !-- if topography is read from ASCII file as no distinction between buildings
@@ -1038,6 +1045,6 @@
 
 !
-!--       Finally, determine maximumum terrain height occupied by the 
-!--       respective building. 
+!--       Determine maximumum terrain height occupied by the respective 
+!--       building and temporalily store on oro_max 
           ALLOCATE( oro_max_l(1:SIZE(build_ids_final)) )
           ALLOCATE( oro_max(1:SIZE(build_ids_final))   )
@@ -1059,9 +1066,18 @@
           oro_max = oro_max_l
 #endif
+!
+!--       Finally, determine discrete grid height of maximum orography occupied
+!--       by a building. Use all-or-nothing approach, i.e. a grid box is either
+          oro_max_l = 0.0
+          DO  nr = 1, SIZE(build_ids_final)
+             DO  k = nzb, nzt
+                IF ( zu(k) - ocean_offset <= oro_max(nr) )                     &
+                   oro_max_l = zw(k) - ocean_offset
+             ENDDO
+             oro_max = oro_max_l
+          ENDDO
        ENDIF
 !
 !--    Map orography as well as buildings onto grid. 
-!--    In case of ocean simulations, add an offset.  
-       ocean_offset = MERGE( zw(0), 0.0_wp, ocean )
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -1080,5 +1096,5 @@
                     topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
                     topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
-                    topo_top_index = topo_top_index + 1
+                    topo_top_index = k ! topo_top_index + 1
                 ENDIF
 !
@@ -1116,16 +1132,25 @@
                                      building_id_f%var(j,i) ), DIM = 1 )
 !
-!--                Extend building down to the terrain surface.
-                   k2 = topo_top_index
+!--                Extend building down to the terrain surface, i.e. fill-up
+!--                surface irregularities below a building. Note, oro_max
+!--                is already a discrete height according to the all-or-nothing
+!--                approach, i.e. grid box is either topography or atmosphere, 
+!--                terrain top is defined at upper bound of the grid box.
+!--                Hence, check for zw in this case. 
                    DO k = topo_top_index + 1, nzt + 1      
-                      IF ( zu(k) - ocean_offset <= oro_max(nr) )  THEN
+                      IF ( zw(k) - ocean_offset <= oro_max(nr) )  THEN
                          topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
                          topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
-                         k2             = k2 + 1
                       ENDIF
-                   ENDDO    
-                   topo_top_index = k2       
-!
-!--                Now, map building on top. 
+                   ENDDO        
+!
+!--                After surface irregularities are smoothen, determine lower
+!--                start index where building starts. 
+                   DO  k = nzb, nzt
+                      IF ( zw(k) - ocean_offset <= oro_max(nr) )               &
+                         topo_top_index = k
+                   ENDDO
+!
+!--                Finally, map building on top.
                    k2 = 0
                    DO k = topo_top_index, nzt + 1
Index: /palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 3050)
+++ /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 3051)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfix in surface-element loops for pavement surfaces 
+! 
+! 3045 2018-05-28 07:55:41Z Giersch
 ! Error messages revised
 ! 
@@ -3788,6 +3791,6 @@
                    surf_lsm_v(l)%emissivity(ind_pav_green,m)  =                &
                                               pavement_pars(ind_p_emis,st)    
-
-                   DO  k = nzb_soil, surf_lsm_h%nzt_pavement(m)
+                                              
+                   DO  k = nzb_soil, surf_lsm_v(l)%nzt_pavement(m)
                       surf_lsm_v(l)%lambda_h_def(k,m)    =                     &
                                     pavement_subsurface_pars_1(k,pavement_type)                       
@@ -4410,5 +4413,5 @@
                 DO  k = nzb_soil, nzt_soil
                    IF ( surf_lsm_v(l)%pavement_surface(m)  .AND.               &
-                        k <= surf_lsm_h%nzt_pavement(m) )  THEN
+                        k <= surf_lsm_v(l)%nzt_pavement(m) )  THEN
                       surf_lsm_v(l)%root_fr(k,m) = 0.0_wp
                    ELSE
@@ -4692,7 +4695,7 @@
           ALLOCATE ( surf_lsm_v(l)%water_surface(1:surf_lsm_v(l)%ns)       )
 
-          surf_lsm_v(l)%water_surface     = .FALSE.
-          surf_lsm_v(l)%pavement_surface  = .FALSE.
-          surf_lsm_v(l)%vegetation_surface   = .FALSE. 
+          surf_lsm_v(l)%water_surface       = .FALSE.
+          surf_lsm_v(l)%pavement_surface    = .FALSE.
+          surf_lsm_v(l)%vegetation_surface  = .FALSE. 
           
 
Index: /palm/trunk/SOURCE/netcdf_data_input_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3050)
+++ /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3051)
@@ -25,8 +25,15 @@
 ! -----------------
 ! $Id$
-! Error messages revised
+! - Speed-up NetCDF input
+! - Revise input routines and remove NetCDF input via IO-blocks since this is 
+!   not working in parallel mode in case blocking collective read operations 
+!   are done 
+! - Temporarily revoke renaming of input variables in dynamic driver (tend_ug, 
+!   tend_vg, zsoil) in order to keep dynamic input file working with current 
+!   model version
+! - More detailed error messages created
 ! 
 ! 3045 2018-05-28 07:55:41Z Giersch
-! Error message revised
+! Error messages revised
 ! 
 ! 3041 2018-05-25 10:39:54Z gronemeier
@@ -563,5 +570,6 @@
        MODULE PROCEDURE get_variable_3d_int8
        MODULE PROCEDURE get_variable_3d_real
-       MODULE PROCEDURE get_variable_3d_real_v
+       MODULE PROCEDURE get_variable_3d_real_dynamic
+!        MODULE PROCEDURE get_variable_3d_real_v
        MODULE PROCEDURE get_variable_4d_real
     END INTERFACE get_variable
@@ -647,47 +655,46 @@
        IF ( .NOT. input_pids_static )  RETURN
 
-       DO  ii = 0, io_blocks-1
-          IF ( ii == io_group )  THEN
 #if defined ( __netcdf )
 !
-!--          Open file in read-only mode
-             CALL open_read_file( TRIM( input_file_static ) //                 &
-                                  TRIM( coupling_char ), id_mod )
-!
-!--          Read global attributes
-             CALL get_attribute( id_mod, input_file_atts%origin_lat_char,      &
-                                 input_file_atts%origin_lat, .TRUE. )
-
-             CALL get_attribute( id_mod, input_file_atts%origin_lon_char,      &
-                                 input_file_atts%origin_lon, .TRUE. )
-
-             CALL get_attribute( id_mod, input_file_atts%origin_time_char,     &
-                                 input_file_atts%origin_time, .TRUE. )
-
-             CALL get_attribute( id_mod, input_file_atts%origin_x_char,        &
-                                 input_file_atts%origin_x, .TRUE. )
-
-             CALL get_attribute( id_mod, input_file_atts%origin_y_char,        &
-                                 input_file_atts%origin_y, .TRUE. )
-
-             CALL get_attribute( id_mod, input_file_atts%origin_z_char,        &
-                                 input_file_atts%origin_z, .TRUE. )
-
-             CALL get_attribute( id_mod, input_file_atts%rotation_angle_char,  &
-                                 input_file_atts%rotation_angle, .TRUE. )
-
-!
-!--          Finally, close input file
-             CALL close_input_file( id_mod )
-!
-!--          Copy lon/lat values
-             init_model%latitude = input_file_atts%origin_lat
-             init_model%longitude = input_file_atts%origin_lon
+!--    Open file in read-only mode
+       CALL open_read_file( TRIM( input_file_static ) //                       &
+                            TRIM( coupling_char ), id_mod )
+!
+!--    Read global attributes
+       CALL get_attribute( id_mod, input_file_atts%origin_lat_char,            &
+                           input_file_atts%origin_lat, .TRUE. )
+
+       CALL get_attribute( id_mod, input_file_atts%origin_lon_char,            &
+                           input_file_atts%origin_lon, .TRUE. )
+
+       CALL get_attribute( id_mod, input_file_atts%origin_time_char,           &
+                           input_file_atts%origin_time, .TRUE. )
+
+       CALL get_attribute( id_mod, input_file_atts%origin_x_char,              &
+                           input_file_atts%origin_x, .TRUE. )
+
+       CALL get_attribute( id_mod, input_file_atts%origin_y_char,              &
+                           input_file_atts%origin_y, .TRUE. )
+
+       CALL get_attribute( id_mod, input_file_atts%origin_z_char,              &
+                           input_file_atts%origin_z, .TRUE. )
+
+       CALL get_attribute( id_mod, input_file_atts%rotation_angle_char,        &
+                           input_file_atts%rotation_angle, .TRUE. )
+
+!
+!--    Finally, close input file
+       CALL close_input_file( id_mod )
 #endif
-          ENDIF
-#if defined( __parallel )
-          CALL MPI_BARRIER( comm2d, ierr )
-#endif
-       ENDDO
+!
+!--    Copy latitude, longitude, origin_z, rotation angle on init type
+       init_model%latitude        = input_file_atts%origin_lat
+       init_model%longitude       = input_file_atts%origin_lon
+       init_model%origin_time     = input_file_atts%origin_time  
+       init_model%origin_x        = input_file_atts%origin_x
+       init_model%origin_y        = input_file_atts%origin_y
+       init_model%origin_z        = input_file_atts%origin_z  
+       init_model%rotation_angle  = input_file_atts%rotation_angle  
+            
 !
 !--    In case of nested runs, each model domain might have different longitude
@@ -697,5 +704,5 @@
 !--    synchronization is required already here.
 #if defined( __parallel )
-       CALL MPI_BCAST( init_model%latitude, 1, MPI_REAL, 0,                    &
+       CALL MPI_BCAST( init_model%latitude,  1, MPI_REAL, 0,                   &
                        MPI_COMM_WORLD, ierr )
        CALL MPI_BCAST( init_model%longitude, 1, MPI_REAL, 0,                   &
@@ -725,8 +732,5 @@
        CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names  !< variable names in static input file
 
-       INTEGER(iwp) ::  i         !< running index along x-direction
-       INTEGER(iwp) ::  ii        !< running index for IO blocks
        INTEGER(iwp) ::  id_surf   !< NetCDF id of input file
-       INTEGER(iwp) ::  j         !< running index along y-direction
        INTEGER(iwp) ::  k         !< running index along z-direction
        INTEGER(iwp) ::  k2        !< running index
@@ -751,109 +755,99 @@
 !--    Read plant canopy variables.
        IF ( plant_canopy )  THEN
-          DO  ii = 0, io_blocks-1
-             IF ( ii == io_group )  THEN
 #if defined ( __netcdf )
 !
-!--             Open file in read-only mode
-                CALL open_read_file( TRIM( input_file_static ) //              &
-                                     TRIM( coupling_char ) , id_surf )
-!
-!--             At first, inquire all variable names.
-!--             This will be used to check whether an optional input variable
-!--             exist or not.
-                CALL inquire_num_variables( id_surf, num_vars )
-
-                ALLOCATE( var_names(1:num_vars) )
-                CALL inquire_variable_names( id_surf, var_names )
-
-!
-!--             Read leaf area density - resolved vegetation
-                IF ( check_existence( var_names, 'lad' ) )  THEN
-                   leaf_area_density_f%from_file = .TRUE.
-                   CALL get_attribute( id_surf, char_fill,                     &
-                                       leaf_area_density_f%fill,               &
-                                       .FALSE., 'lad' )
-!
-!--                Inquire number of vertical vegetation layer
-                   CALL get_dimension_length( id_surf, leaf_area_density_f%nz, &
-                                              'zlad' )
-!
-!--                Allocate variable for leaf-area density
-                   ALLOCATE( leaf_area_density_f%var(                          &
-                                                   0:leaf_area_density_f%nz-1, &
-                                                   nys:nyn,nxl:nxr) )
-
-                   CALL get_variable( id_surf, 'lad',      &
-                                      nxl, nxr, nys, nyn,                &
-                                      leaf_area_density_f%var(:,nys:nyn, nxl:nxr) )
-
-                ELSE
-                   leaf_area_density_f%from_file = .FALSE.
-                ENDIF
-
-!
-!--             Read basal area density - resolved vegetation
-                IF ( check_existence( var_names, 'bad' ) )  THEN
-                   basal_area_density_f%from_file = .TRUE.
-                   CALL get_attribute( id_surf, char_fill,                     &
-                                       basal_area_density_f%fill,              &
-                                       .FALSE., 'bad' )
-!
-!--                Inquire number of vertical vegetation layer
-                   CALL get_dimension_length( id_surf,                         &
-                                              basal_area_density_f%nz,         &
-                                              'zlad' )
-!
-!--                Allocate variable
-                   ALLOCATE( basal_area_density_f%var(                         &
-                                                  0:basal_area_density_f%nz-1, &
-                                                  nys:nyn,nxl:nxr) )
-
-                   CALL get_variable( id_surf, 'bad',                    &
-                                      nxl, nxr, nys, nyn,                &
-                                      basal_area_density_f%var(:,nys:nyn, nxl:nxr) )
-                ELSE
-                   basal_area_density_f%from_file = .FALSE.
-                ENDIF
-
-!
-!--             Read root area density - resolved vegetation
-                IF ( check_existence( var_names, 'root_area_dens_r' ) )  THEN
-                   root_area_density_lad_f%from_file = .TRUE.
-                   CALL get_attribute( id_surf, char_fill,                     &
-                                       root_area_density_lad_f%fill,           &
-                                       .FALSE., 'root_area_dens_r' )
-!
-!--                Inquire number of vertical soil layers
-                   CALL get_dimension_length( id_surf,                         &
-                                              root_area_density_lad_f%nz,      &
-                                              'zsoil' )
-!
-!--                Allocate variable
-                   ALLOCATE( root_area_density_lad_f%var                       &
-                                               (0:root_area_density_lad_f%nz-1,&
+!--       Open file in read-only mode
+          CALL open_read_file( TRIM( input_file_static ) //                    &
+                               TRIM( coupling_char ) , id_surf )
+!
+!--       At first, inquire all variable names.
+!--       This will be used to check whether an optional input variable
+!--       exist or not.
+          CALL inquire_num_variables( id_surf, num_vars )
+
+          ALLOCATE( var_names(1:num_vars) )
+          CALL inquire_variable_names( id_surf, var_names )
+
+!
+!--       Read leaf area density - resolved vegetation
+          IF ( check_existence( var_names, 'lad' ) )  THEN
+             leaf_area_density_f%from_file = .TRUE.
+             CALL get_attribute( id_surf, char_fill,                           &
+                                 leaf_area_density_f%fill,                     &
+                                 .FALSE., 'lad' )
+!
+!--          Inquire number of vertical vegetation layer
+             CALL get_dimension_length( id_surf, leaf_area_density_f%nz,       &
+                                        'zlad' )
+!
+!--          Allocate variable for leaf-area density
+             ALLOCATE( leaf_area_density_f%var( 0:leaf_area_density_f%nz-1,    &
                                                 nys:nyn,nxl:nxr) )
 
-                   CALL get_variable( id_surf, 'root_area_dens_r',  &
-                                      nxl, nxr, nys, nyn,           &
-                                      root_area_density_lad_f%var(:,nys:nyn, nxl:nxr) )
-                ELSE
-                   root_area_density_lad_f%from_file = .FALSE.
-                ENDIF
-!
-!--             Finally, close input file
-                CALL close_input_file( id_surf )
+             CALL get_variable( id_surf, 'lad', leaf_area_density_f%var,       &
+                                nxl, nxr, nys, nyn,                            &
+                                0, leaf_area_density_f%nz-1 )
+
+          ELSE
+             leaf_area_density_f%from_file = .FALSE.
+          ENDIF
+
+!
+!--       Read basal area density - resolved vegetation
+          IF ( check_existence( var_names, 'bad' ) )  THEN
+             basal_area_density_f%from_file = .TRUE.
+             CALL get_attribute( id_surf, char_fill,                           &
+                                 basal_area_density_f%fill,                    &
+                                 .FALSE., 'bad' )
+!
+!--          Inquire number of vertical vegetation layer
+             CALL get_dimension_length( id_surf, basal_area_density_f%nz,      &
+                                        'zlad' )
+!
+!--          Allocate variable
+             ALLOCATE( basal_area_density_f%var(0:basal_area_density_f%nz-1,   &
+                                                nys:nyn,nxl:nxr) )
+
+             CALL get_variable( id_surf, 'bad', basal_area_density_f%var,      &
+                                nxl, nxr, nys, nyn,                            &
+                                0,  basal_area_density_f%nz-1 )
+          ELSE
+             basal_area_density_f%from_file = .FALSE.
+          ENDIF
+
+!
+!--       Read root area density - resolved vegetation
+          IF ( check_existence( var_names, 'root_area_dens_r' ) )  THEN
+             root_area_density_lad_f%from_file = .TRUE.
+             CALL get_attribute( id_surf, char_fill,                           &
+                                 root_area_density_lad_f%fill,                 &
+                                 .FALSE., 'root_area_dens_r' )
+!
+!--          Inquire number of vertical soil layers
+             CALL get_dimension_length( id_surf,                               &
+                                        root_area_density_lad_f%nz,            &
+                                        'zsoil' )
+!
+!--          Allocate variable
+             ALLOCATE( root_area_density_lad_f%var                             &
+                                         (0:root_area_density_lad_f%nz-1,      &
+                                          nys:nyn,nxl:nxr) )
+
+             CALL get_variable( id_surf, 'root_area_dens_r',                   &
+                                root_area_density_lad_f%var,                   &
+                                nxl, nxr, nys, nyn,                            &
+                                0,  root_area_density_lad_f%nz-1 )
+          ELSE
+             root_area_density_lad_f%from_file = .FALSE.
+          ENDIF
+!
+!--       Finally, close input file
+          CALL close_input_file( id_surf )
 #endif
-             ENDIF
-#if defined( __parallel )
-             CALL MPI_BARRIER( comm2d, ierr )
-#endif
-          ENDDO
-!
-!--       Deallocate variable list. Will be re-allocated in case further
-!--       variables are read from file.
-          IF ( ALLOCATED( var_names ) )  DEALLOCATE( var_names )
-
-       ENDIF
+       ENDIF
+!
+!--    Deallocate variable list. Will be re-allocated in case further
+!--    variables are read from file.
+       IF ( ALLOCATED( var_names ) )  DEALLOCATE( var_names )
 !
 !--    Skip the following if no land-surface or urban-surface module are
@@ -865,487 +859,439 @@
        var_exchange_real = 0.0_wp
 
-       DO  ii = 0, io_blocks-1
-          IF ( ii == io_group )  THEN
 #if defined ( __netcdf )
 !
-!--          Open file in read-only mode
-             CALL open_read_file( TRIM( input_file_static ) //                 &
-                                  TRIM( coupling_char ) , id_surf )
-
-!
-!--          Inquire all variable names.
-!--          This will be used to check whether an optional input variable exist
-!--          or not.
-             CALL inquire_num_variables( id_surf, num_vars )
-
-             ALLOCATE( var_names(1:num_vars) )
-             CALL inquire_variable_names( id_surf, var_names )
-
-!
-!--          Read vegetation type and required attributes
-             IF ( check_existence( var_names, 'vegetation_type' ) )  THEN
-                vegetation_type_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill,                        &
-                                    vegetation_type_f%fill,                    &
-                                    .FALSE., 'vegetation_type' )
-!
-!--             PE-wise reading of 2D vegetation type.
-                ALLOCATE ( vegetation_type_f%var(nys:nyn,nxl:nxr)  )
-
-                DO  i = nxl, nxr
-                   CALL get_variable( id_surf, 'vegetation_type',              &
-                                      i, vegetation_type_f%var(:,i) )
-                ENDDO
-             ELSE
-                vegetation_type_f%from_file = .FALSE.
-             ENDIF
-
-!
-!--          Read soil type and required attributes
-             IF ( check_existence( var_names, 'soil_type' ) )  THEN
-                   soil_type_f%from_file = .TRUE.
-!
-!--             Note, lod is currently not on file; skip for the moment
-!                 CALL get_attribute( id_surf, char_lod,                       &
-!                                            soil_type_f%lod,                  &
-!                                            .FALSE., 'soil_type' )
-                CALL get_attribute( id_surf, char_fill,                        &
-                                    soil_type_f%fill,                          &
-                                    .FALSE., 'soil_type' )
-
-                IF ( soil_type_f%lod == 1 )  THEN
-!
-!--                PE-wise reading of 2D soil type.
-                   ALLOCATE ( soil_type_f%var_2d(nys:nyn,nxl:nxr)  )
-                   DO  i = nxl, nxr
-                      CALL get_variable( id_surf, 'soil_type',                 &
-                                         i, soil_type_f%var_2d(:,i) )
-                   ENDDO
-                ELSEIF ( soil_type_f%lod == 2 )  THEN
-!
-!--                Obtain number of soil layers from file.
-                   CALL get_dimension_length( id_surf, nz_soil, 'zsoil' )
-!
-!--                PE-wise reading of 3D soil type.
-                   ALLOCATE ( soil_type_f%var_3d(0:nz_soil,nys:nyn,nxl:nxr) )
-                   DO  i = nxl, nxr
-                      DO  j = nys, nyn
-                         CALL get_variable( id_surf, 'soil_type', i, j,        &
-                                            soil_type_f%var_3d(:,j,i) )
-                      ENDDO
-                   ENDDO
-                ENDIF
-             ELSE
-                soil_type_f%from_file = .FALSE.
-             ENDIF
-
-!
-!--          Read pavement type and required attributes
-             IF ( check_existence( var_names, 'pavement_type' ) )  THEN
-                pavement_type_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill,                        &
-                                    pavement_type_f%fill, .FALSE.,             &
-                                    'pavement_type' )
-!
-!--             PE-wise reading of 2D pavement type.
-                ALLOCATE ( pavement_type_f%var(nys:nyn,nxl:nxr)  )
-                DO  i = nxl, nxr
-                   CALL get_variable( id_surf, 'pavement_type',                &
-                                      i, pavement_type_f%var(:,i) )
-                ENDDO
-             ELSE
-                pavement_type_f%from_file = .FALSE.
-             ENDIF
-
-!
-!--          Read water type and required attributes
-             IF ( check_existence( var_names, 'water_type' ) )  THEN
-                water_type_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill, water_type_f%fill,     &
-                                    .FALSE., 'water_type' )
-!
-!--             PE-wise reading of 2D water type.
-                ALLOCATE ( water_type_f%var(nys:nyn,nxl:nxr)  )
-                DO  i = nxl, nxr
-                   CALL get_variable( id_surf, 'water_type', i,                &
-                                      water_type_f%var(:,i) )
-                ENDDO
-             ELSE
-                water_type_f%from_file = .FALSE.
-             ENDIF
-!
-!--          Read surface fractions and related information
-             IF ( check_existence( var_names, 'surface_fraction' ) )  THEN
-                surface_fraction_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill,                        &
-                                    surface_fraction_f%fill,                   &
-                                    .FALSE., 'surface_fraction' )
-!
-!--             Inquire number of surface fractions
-                CALL get_dimension_length( id_surf,                            &
-                                           surface_fraction_f%nf,              &
-                                           'nsurface_fraction' )
-!
-!--             Allocate dimension array and input array for surface fractions
-                ALLOCATE( surface_fraction_f%nfracs(0:surface_fraction_f%nf-1) )
-                ALLOCATE( surface_fraction_f%frac(0:surface_fraction_f%nf-1,   &
-                                                  nys:nyn,nxl:nxr) )
-!
-!--             Get dimension of surface fractions
-                CALL get_variable( id_surf, 'nsurface_fraction',               &
-                                   surface_fraction_f%nfracs )
-!
-!--             Read surface fractions
-                CALL get_variable( id_surf, 'surface_fraction',     &
-                                   nxl, nxr, nys, nyn,              &
-                                   surface_fraction_f%frac(:,nys:nyn, nxl:nxr))
-             ELSE
-                surface_fraction_f%from_file = .FALSE.
-             ENDIF
-!
-!--          Read building parameters and related information
-             IF ( check_existence( var_names, 'building_pars' ) )  THEN
-                building_pars_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill,                        &
-                                    building_pars_f%fill,                      &
-                                    .FALSE., 'building_pars' )
-!
-!--             Inquire number of building parameters
-                CALL get_dimension_length( id_surf,                            &
-                                           building_pars_f%np,                 &
-                                           'nbuilding_pars' )
-!
-!--             Allocate dimension array and input array for building parameters
-                ALLOCATE( building_pars_f%pars(0:building_pars_f%np-1) )
-                ALLOCATE( building_pars_f%pars_xy(0:building_pars_f%np-1,      &
-                                                  nys:nyn,nxl:nxr) )
-!
-!--             Get dimension of building parameters
-                CALL get_variable( id_surf, 'nbuilding_pars',                  &
-                                   building_pars_f%pars )
-!
-!--             Read building_pars
-                CALL get_variable( id_surf, 'building_pars',                   &
-                                   nxl, nxr, nys, nyn,                         &
-                                   building_pars_f%pars_xy(:,nys:nyn, nxl:nxr) )
-             ELSE
-                building_pars_f%from_file = .FALSE.
-             ENDIF
-
-!
-!--          Read albedo type and required attributes
-             IF ( check_existence( var_names, 'albedo_type' ) )  THEN
-                albedo_type_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill, albedo_type_f%fill,    &
-                                    .FALSE.,  'albedo_type' )
-!
-!--             PE-wise reading of 2D water type.
-                ALLOCATE ( albedo_type_f%var(nys:nyn,nxl:nxr)  )
-                DO  i = nxl, nxr
-                   CALL get_variable( id_surf, 'albedo_type',                  &
-                                      i, albedo_type_f%var(:,i) )
-                ENDDO
-             ELSE
-                albedo_type_f%from_file = .FALSE.
-             ENDIF
-!
-!--          Read albedo parameters and related information
-             IF ( check_existence( var_names, 'albedo_pars' ) )  THEN
-                albedo_pars_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill, albedo_pars_f%fill,    &
-                                    .FALSE., 'albedo_pars' )
-!
-!--             Inquire number of albedo parameters
-                CALL get_dimension_length( id_surf, albedo_pars_f%np,          &
-                                           'nalbedo_pars' )
-!
-!--             Allocate dimension array and input array for albedo parameters
-                ALLOCATE( albedo_pars_f%pars(0:albedo_pars_f%np-1) )
-                ALLOCATE( albedo_pars_f%pars_xy(0:albedo_pars_f%np-1,          &
-                                                nys:nyn,nxl:nxr) )
-!
-!--             Get dimension of albedo parameters
-                CALL get_variable( id_surf, 'nalbedo_pars', albedo_pars_f%pars )
-
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      CALL get_variable( id_surf, 'albedo_pars', i, j,         &
-                                         albedo_pars_f%pars_xy(:,j,i) )
-                   ENDDO
-                ENDDO
-             ELSE
-                albedo_pars_f%from_file = .FALSE.
-             ENDIF
-
-!
-!--          Read pavement parameters and related information
-             IF ( check_existence( var_names, 'pavement_pars' ) )  THEN
-                pavement_pars_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill,                        &
-                                    pavement_pars_f%fill,                      &
-                                    .FALSE., 'pavement_pars' )
-!
-!--             Inquire number of pavement parameters
-                CALL get_dimension_length( id_surf, pavement_pars_f%np,        &
-                                           'npavement_pars' )
-!
-!--             Allocate dimension array and input array for pavement parameters
-                ALLOCATE( pavement_pars_f%pars(0:pavement_pars_f%np-1) )
-                ALLOCATE( pavement_pars_f%pars_xy(0:pavement_pars_f%np-1,      &
-                                                  nys:nyn,nxl:nxr) )
-!
-!--             Get dimension of pavement parameters
-                CALL get_variable( id_surf, 'npavement_pars',                  &
-                                   pavement_pars_f%pars )
-
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      CALL get_variable( id_surf, 'pavement_pars', i, j,       &
-                                         pavement_pars_f%pars_xy(:,j,i) )
-                   ENDDO
-                ENDDO
-             ELSE
-                pavement_pars_f%from_file = .FALSE.
-             ENDIF
-
-!
-!--          Read pavement subsurface parameters and related information
-             IF ( check_existence( var_names, 'pavement_subsurface_pars' ) )   &
-             THEN
-                pavement_subsurface_pars_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill,                        &
-                                    pavement_subsurface_pars_f%fill,           &
-                                    .FALSE., 'pavement_subsurface_pars' )
-!
-!--             Inquire number of parameters
-                CALL get_dimension_length( id_surf,                            &
-                                           pavement_subsurface_pars_f%np,      &
-                                           'npavement_subsurface_pars' )
-!
-!--             Inquire number of soil layers
-                CALL get_dimension_length( id_surf,                            &
-                                           pavement_subsurface_pars_f%nz,      &
-                                           'zsoil' )
-!
-!--             Allocate dimension array and input array for pavement parameters
-                ALLOCATE( pavement_subsurface_pars_f%pars                      &
-                                  (0:pavement_subsurface_pars_f%np-1) )
-                ALLOCATE( pavement_subsurface_pars_f%pars_xyz                  &
-                                  (0:pavement_subsurface_pars_f%np-1,          &
-                                   0:pavement_subsurface_pars_f%nz-1,          &
-                                   nys:nyn,nxl:nxr) )
-!
-!--             Get dimension of pavement parameters
-                CALL get_variable( id_surf, 'npavement_subsurface_pars',       &
-                                   pavement_subsurface_pars_f%pars )
-
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      CALL get_variable(                                       &
-                                  id_surf, 'pavement_subsurface_pars',         &
-                                  i, j,                                        &
-                                  pavement_subsurface_pars_f%pars_xyz(:,:,j,i),&
-                                  pavement_subsurface_pars_f%nz,               &
-                                  pavement_subsurface_pars_f%np )
-                   ENDDO
-                ENDDO
-             ELSE
-                pavement_subsurface_pars_f%from_file = .FALSE.
-             ENDIF
-
-
-!
-!--          Read vegetation parameters and related information
-             IF ( check_existence( var_names, 'vegetation_pars' ) )  THEN
-                vegetation_pars_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill,                        &
-                                    vegetation_pars_f%fill,                    &
-                                    .FALSE.,  'vegetation_pars' )
-!
-!--             Inquire number of vegetation parameters
-                CALL get_dimension_length( id_surf, vegetation_pars_f%np,      &
-                                           'nvegetation_pars' )
-!
-!--             Allocate dimension array and input array for surface fractions
-                ALLOCATE( vegetation_pars_f%pars(0:vegetation_pars_f%np-1) )
-                ALLOCATE( vegetation_pars_f%pars_xy(0:vegetation_pars_f%np-1,  &
-                                                    nys:nyn,nxl:nxr) )
-!
-!--             Get dimension of the parameters
-                CALL get_variable( id_surf, 'nvegetation_pars',                &
-                                   vegetation_pars_f%pars )
-
-                CALL get_variable( id_surf, 'vegetation_pars',                 &
-                                   nxl, nxr, nys, nyn,                         &
-                                   vegetation_pars_f%pars_xy(:,nys:nyn,nxl:nxr ) )
-             ELSE
-                vegetation_pars_f%from_file = .FALSE.
-             ENDIF
-
-!
-!--          Read root parameters/distribution and related information
-             IF ( check_existence( var_names, 'soil_pars' ) )  THEN
-                soil_pars_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill,                        &
-                                    soil_pars_f%fill,                          &
-                                    .FALSE., 'soil_pars' )
-
-                CALL get_attribute( id_surf, char_lod,                         &
-                                    soil_pars_f%lod,                           &
-                                    .FALSE., 'soil_pars' )
-
-!
-!--             Inquire number of soil parameters
-                CALL get_dimension_length( id_surf,                            &
-                                           soil_pars_f%np,                     &
-                                           'nsoil_pars' )
-!
-!--             Read parameters array
-                ALLOCATE( soil_pars_f%pars(0:soil_pars_f%np-1) )
-                CALL get_variable( id_surf, 'nsoil_pars', soil_pars_f%pars )
-
-!
-!--             In case of level of detail 2, also inquire number of vertical
-!--             soil layers, allocate memory and read the respective dimension
-                IF ( soil_pars_f%lod == 2 )  THEN
-                   CALL get_dimension_length( id_surf, soil_pars_f%nz, 'zsoil' )
-
-                   ALLOCATE( soil_pars_f%layers(0:soil_pars_f%nz-1) )
-                   CALL get_variable( id_surf, 'zsoil', soil_pars_f%layers )
-
-                ENDIF
-
-!
-!--             Read soil parameters, depending on level of detail
-                IF ( soil_pars_f%lod == 1 )  THEN
-                   ALLOCATE( soil_pars_f%pars_xy(0:soil_pars_f%np-1,           &
-                                                 nys:nyn,nxl:nxr) )
-                   DO  i = nxl, nxr
-                      DO  j = nys, nyn
-                         CALL get_variable( id_surf, 'soil_pars', i, j,        &
-                                            soil_pars_f%pars_xy(:,j,i) )
-                      ENDDO
-                   ENDDO
-
-                ELSEIF ( soil_pars_f%lod == 2 )  THEN
-                   ALLOCATE( soil_pars_f%pars_xyz(0:soil_pars_f%np-1,          &
-                                                  0:soil_pars_f%nz-1,          &
-                                                  nys:nyn,nxl:nxr) )
-                   DO  i = nxl, nxr
-                      DO  j = nys, nyn
-                         CALL get_variable( id_surf, 'soil_pars', i, j,        &
-                                            soil_pars_f%pars_xyz(:,:,j,i),     &
-                                            soil_pars_f%nz, soil_pars_f%np )
-                      ENDDO
-                   ENDDO
-                ENDIF
-             ELSE
-                soil_pars_f%from_file = .FALSE.
-             ENDIF
-
-!
-!--          Read water parameters and related information
-             IF ( check_existence( var_names, 'water_pars' ) )  THEN
-                water_pars_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill,                        &
-                                    water_pars_f%fill,                         &
-                                    .FALSE., 'water_pars' )
-!
-!--             Inquire number of water parameters
-                CALL get_dimension_length( id_surf,                            &
-                                           water_pars_f%np,                    &
-                                           'nwater_pars' )
-!
-!--             Allocate dimension array and input array for water parameters
-                ALLOCATE( water_pars_f%pars(0:water_pars_f%np-1) )
-                ALLOCATE( water_pars_f%pars_xy(0:water_pars_f%np-1,            &
-                                               nys:nyn,nxl:nxr) )
-!
-!--             Get dimension of water parameters
-                CALL get_variable( id_surf, 'nwater_pars', water_pars_f%pars )
-
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      CALL get_variable( id_surf, 'water_pars', i, j,          &
-                                         water_pars_f%pars_xy(:,j,i) )
-                   ENDDO
-                ENDDO
-             ELSE
-                water_pars_f%from_file = .FALSE.
-             ENDIF
-!
-!--          Read root area density - parametrized vegetation
-             IF ( check_existence( var_names, 'root_area_dens_s' ) )  THEN
-                root_area_density_lsm_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill,                        &
-                                    root_area_density_lsm_f%fill,              &
-                                    .FALSE., 'root_area_dens_s' )
-!
-!--             Obtain number of soil layers from file and allocate variable
-                CALL get_dimension_length( id_surf, root_area_density_lsm_f%nz,&
-                                           'zsoil' )
-                ALLOCATE( root_area_density_lsm_f%var                          &
-                                              (0:root_area_density_lsm_f%nz-1, &
-                                               nys:nyn,nxl:nxr) )
-
-!
-!--             Read root-area density
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      CALL get_variable( id_surf, 'root_area_dens_s',          &
-                                         i, j,                                 &
-                                         root_area_density_lsm_f%var(:,j,i) )
-                   ENDDO
-                ENDDO
-
-             ELSE
-                root_area_density_lsm_f%from_file = .FALSE.
-             ENDIF
-!
-!--          Read street type and street crossing
-             IF ( check_existence( var_names, 'street_type' ) )  THEN
-                street_type_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill,                        &
-                                    street_type_f%fill, .FALSE.,               &
-                                    'street_type' )
-!
-!--             PE-wise reading of 2D pavement type.
-                ALLOCATE ( street_type_f%var(nys:nyn,nxl:nxr)  )
-                DO  i = nxl, nxr
-                   CALL get_variable( id_surf, 'street_type',                  &
-                                      i, street_type_f%var(:,i) )
-                ENDDO
-             ELSE
-                street_type_f%from_file = .FALSE.
-             ENDIF
-
-             IF ( check_existence( var_names, 'street_crossing' ) )  THEN
-                street_crossing_f%from_file = .TRUE.
-                CALL get_attribute( id_surf, char_fill,                        &
-                                    street_crossing_f%fill, .FALSE.,           &
-                                    'street_crossing' )
-!
-!--             PE-wise reading of 2D pavement type.
-                ALLOCATE ( street_crossing_f%var(nys:nyn,nxl:nxr)  )
-                DO  i = nxl, nxr
-                   CALL get_variable( id_surf, 'street_crossing',              &
-                                      i, street_crossing_f%var(:,i) )
-                ENDDO
-             ELSE
-                street_crossing_f%from_file = .FALSE.
-             ENDIF
-!
-!--          Still missing: root_resolved and building_surface_pars.
-!--          Will be implemented as soon as they are available.
-
-!
-!--          Finally, close input file
-             CALL close_input_file( id_surf )
+!--    Open file in read-only mode
+       CALL open_read_file( TRIM( input_file_static ) //                       &
+                            TRIM( coupling_char ) , id_surf )
+!
+!--    Inquire all variable names.
+!--    This will be used to check whether an optional input variable exist
+!--    or not.
+       CALL inquire_num_variables( id_surf, num_vars )
+
+       ALLOCATE( var_names(1:num_vars) )
+       CALL inquire_variable_names( id_surf, var_names )
+!
+!--    Read vegetation type and required attributes
+       IF ( check_existence( var_names, 'vegetation_type' ) )  THEN
+          vegetation_type_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill,                              &
+                              vegetation_type_f%fill,                          &
+                              .FALSE., 'vegetation_type' )
+
+          ALLOCATE ( vegetation_type_f%var(nys:nyn,nxl:nxr)  )
+
+          CALL get_variable( id_surf, 'vegetation_type',                       &
+                             vegetation_type_f%var, nxl, nxr, nys, nyn )
+       ELSE
+          vegetation_type_f%from_file = .FALSE.
+       ENDIF
+
+!
+!--    Read soil type and required attributes
+       IF ( check_existence( var_names, 'soil_type' ) )  THEN
+             soil_type_f%from_file = .TRUE.
+!
+!--       Note, lod is currently not on file; skip for the moment
+!           CALL get_attribute( id_surf, char_lod,                       &
+!                                      soil_type_f%lod,                  &
+!                                      .FALSE., 'soil_type' )
+          CALL get_attribute( id_surf, char_fill,                              &
+                              soil_type_f%fill,                                &
+                              .FALSE., 'soil_type' )
+
+          IF ( soil_type_f%lod == 1 )  THEN
+
+             ALLOCATE ( soil_type_f%var_2d(nys:nyn,nxl:nxr)  )
+
+             CALL get_variable( id_surf, 'soil_type', soil_type_f%var_2d,      &
+                                nxl, nxr, nys, nyn )
+
+          ELSEIF ( soil_type_f%lod == 2 )  THEN
+!
+!--          Obtain number of soil layers from file.
+             CALL get_dimension_length( id_surf, nz_soil, 'zsoil' )
+
+             ALLOCATE ( soil_type_f%var_3d(0:nz_soil,nys:nyn,nxl:nxr) )
+
+             CALL get_variable( id_surf, 'soil_type', soil_type_f%var_3d,      &
+                                nxl, nxr, nys, nyn, 0, nz_soil )
+ 
+          ENDIF
+       ELSE
+          soil_type_f%from_file = .FALSE.
+       ENDIF
+
+!
+!--    Read pavement type and required attributes
+       IF ( check_existence( var_names, 'pavement_type' ) )  THEN
+          pavement_type_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill,                              &
+                              pavement_type_f%fill, .FALSE.,                   &
+                              'pavement_type' )
+
+          ALLOCATE ( pavement_type_f%var(nys:nyn,nxl:nxr)  )
+
+          CALL get_variable( id_surf, 'pavement_type', pavement_type_f%var,    &
+                             nxl, nxr, nys, nyn )
+       ELSE
+          pavement_type_f%from_file = .FALSE.
+       ENDIF
+
+!
+!--    Read water type and required attributes
+       IF ( check_existence( var_names, 'water_type' ) )  THEN
+          water_type_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill, water_type_f%fill,           &
+                              .FALSE., 'water_type' )
+
+          ALLOCATE ( water_type_f%var(nys:nyn,nxl:nxr)  )
+
+          CALL get_variable( id_surf, 'water_type', water_type_f%var,          &
+                             nxl, nxr, nys, nyn )
+
+       ELSE
+          water_type_f%from_file = .FALSE.
+       ENDIF
+!
+!--    Read surface fractions and related information
+       IF ( check_existence( var_names, 'surface_fraction' ) )  THEN
+          surface_fraction_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill,                              &
+                              surface_fraction_f%fill,                         &
+                              .FALSE., 'surface_fraction' )
+!
+!--       Inquire number of surface fractions
+          CALL get_dimension_length( id_surf,                                  &
+                                     surface_fraction_f%nf,                    &
+                                     'nsurface_fraction' )
+!
+!--       Allocate dimension array and input array for surface fractions
+          ALLOCATE( surface_fraction_f%nfracs(0:surface_fraction_f%nf-1) )
+          ALLOCATE( surface_fraction_f%frac(0:surface_fraction_f%nf-1,         &
+                                            nys:nyn,nxl:nxr) )
+!
+!--       Get dimension of surface fractions
+          CALL get_variable( id_surf, 'nsurface_fraction',                     &
+                             surface_fraction_f%nfracs )
+!
+!--       Read surface fractions
+          CALL get_variable( id_surf, 'surface_fraction',                      &
+                             surface_fraction_f%frac, nxl, nxr, nys, nyn,      &
+                             0, surface_fraction_f%nf-1 )
+       ELSE
+          surface_fraction_f%from_file = .FALSE.
+       ENDIF
+!
+!--    Read building parameters and related information
+       IF ( check_existence( var_names, 'building_pars' ) )  THEN
+          building_pars_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill,                              &
+                              building_pars_f%fill,                            &
+                              .FALSE., 'building_pars' )
+!
+!--       Inquire number of building parameters
+          CALL get_dimension_length( id_surf,                                  &
+                                     building_pars_f%np,                       &
+                                     'nbuilding_pars' )
+!
+!--       Allocate dimension array and input array for building parameters
+          ALLOCATE( building_pars_f%pars(0:building_pars_f%np-1) )
+          ALLOCATE( building_pars_f%pars_xy(0:building_pars_f%np-1,            &
+                                            nys:nyn,nxl:nxr) )
+!
+!--       Get dimension of building parameters
+          CALL get_variable( id_surf, 'nbuilding_pars',                        &
+                             building_pars_f%pars )
+!
+!--       Read building_pars
+          CALL get_variable( id_surf, 'building_pars',                         &
+                             building_pars_f%pars_xy, nxl, nxr, nys, nyn,      &
+                             0, building_pars_f%np-1 )
+       ELSE
+          building_pars_f%from_file = .FALSE.
+       ENDIF
+
+!
+!--    Read albedo type and required attributes
+       IF ( check_existence( var_names, 'albedo_type' ) )  THEN
+          albedo_type_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill, albedo_type_f%fill,          &
+                              .FALSE.,  'albedo_type' )
+
+          ALLOCATE ( albedo_type_f%var(nys:nyn,nxl:nxr)  )
+          
+          CALL get_variable( id_surf, 'albedo_type', albedo_type_f%var,        &
+                             nxl, nxr, nys, nyn )
+       ELSE
+          albedo_type_f%from_file = .FALSE.
+       ENDIF
+!
+!--    Read albedo parameters and related information
+       IF ( check_existence( var_names, 'albedo_pars' ) )  THEN
+          albedo_pars_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill, albedo_pars_f%fill,          &
+                              .FALSE., 'albedo_pars' )
+!
+!--       Inquire number of albedo parameters
+          CALL get_dimension_length( id_surf, albedo_pars_f%np,                &
+                                     'nalbedo_pars' )
+!
+!--       Allocate dimension array and input array for albedo parameters
+          ALLOCATE( albedo_pars_f%pars(0:albedo_pars_f%np-1) )
+          ALLOCATE( albedo_pars_f%pars_xy(0:albedo_pars_f%np-1,                &
+                                          nys:nyn,nxl:nxr) )
+!
+!--       Get dimension of albedo parameters
+          CALL get_variable( id_surf, 'nalbedo_pars', albedo_pars_f%pars )
+
+          CALL get_variable( id_surf, 'albedo_pars', albedo_pars_f%pars_xy,    &
+                             nxl, nxr, nys, nyn,                               &
+                             0, albedo_pars_f%np-1 )
+       ELSE
+          albedo_pars_f%from_file = .FALSE.
+       ENDIF
+
+!
+!--    Read pavement parameters and related information
+       IF ( check_existence( var_names, 'pavement_pars' ) )  THEN
+          pavement_pars_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill,                              &
+                              pavement_pars_f%fill,                            &
+                              .FALSE., 'pavement_pars' )
+!
+!--       Inquire number of pavement parameters
+          CALL get_dimension_length( id_surf, pavement_pars_f%np,              &
+                                     'npavement_pars' )
+!
+!--       Allocate dimension array and input array for pavement parameters
+          ALLOCATE( pavement_pars_f%pars(0:pavement_pars_f%np-1) )
+          ALLOCATE( pavement_pars_f%pars_xy(0:pavement_pars_f%np-1,            &
+                                            nys:nyn,nxl:nxr) )
+!
+!--       Get dimension of pavement parameters
+          CALL get_variable( id_surf, 'npavement_pars', pavement_pars_f%pars )
+
+          CALL get_variable( id_surf, 'pavement_pars', pavement_pars_f%pars_xy,&
+                             nxl, nxr, nys, nyn,                               &
+                             0, pavement_pars_f%np-1 )
+       ELSE
+          pavement_pars_f%from_file = .FALSE.
+       ENDIF
+
+!
+!--    Read pavement subsurface parameters and related information
+       IF ( check_existence( var_names, 'pavement_subsurface_pars' ) )         &
+       THEN
+          pavement_subsurface_pars_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill,                              &
+                              pavement_subsurface_pars_f%fill,                 &
+                              .FALSE., 'pavement_subsurface_pars' )
+!
+!--       Inquire number of parameters
+          CALL get_dimension_length( id_surf,                                  &
+                                     pavement_subsurface_pars_f%np,            &
+                                     'npavement_subsurface_pars' )
+!
+!--       Inquire number of soil layers
+          CALL get_dimension_length( id_surf,                                  &
+                                     pavement_subsurface_pars_f%nz,            &
+                                     'zsoil' )
+!
+!--       Allocate dimension array and input array for pavement parameters
+          ALLOCATE( pavement_subsurface_pars_f%pars                            &
+                            (0:pavement_subsurface_pars_f%np-1) )
+          ALLOCATE( pavement_subsurface_pars_f%pars_xyz                        &
+                            (0:pavement_subsurface_pars_f%np-1,                &
+                             0:pavement_subsurface_pars_f%nz-1,                &
+                             nys:nyn,nxl:nxr) )
+!
+!--       Get dimension of pavement parameters
+          CALL get_variable( id_surf, 'npavement_subsurface_pars',             &
+                             pavement_subsurface_pars_f%pars )
+
+          CALL get_variable( id_surf, 'pavement_subsurface_pars',              &
+                             pavement_subsurface_pars_f%pars_xyz,              &
+                             nxl, nxr, nys, nyn,                               &
+                             0, pavement_subsurface_pars_f%nz-1,               &
+                             0, pavement_subsurface_pars_f%np-1 )
+       ELSE
+          pavement_subsurface_pars_f%from_file = .FALSE.
+       ENDIF
+
+
+!
+!--    Read vegetation parameters and related information
+       IF ( check_existence( var_names, 'vegetation_pars' ) )  THEN
+          vegetation_pars_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill,                              &
+                              vegetation_pars_f%fill,                          &
+                              .FALSE.,  'vegetation_pars' )
+!
+!--       Inquire number of vegetation parameters
+          CALL get_dimension_length( id_surf, vegetation_pars_f%np,            &
+                                     'nvegetation_pars' )
+!
+!--       Allocate dimension array and input array for surface fractions
+          ALLOCATE( vegetation_pars_f%pars(0:vegetation_pars_f%np-1) )
+          ALLOCATE( vegetation_pars_f%pars_xy(0:vegetation_pars_f%np-1,        &
+                                              nys:nyn,nxl:nxr) )
+!
+!--       Get dimension of the parameters
+          CALL get_variable( id_surf, 'nvegetation_pars',                      &
+                             vegetation_pars_f%pars )
+
+          CALL get_variable( id_surf, 'vegetation_pars',                       &
+                             vegetation_pars_f%pars_xy, nxl, nxr, nys, nyn,    &
+                             0, vegetation_pars_f%np-1 )
+       ELSE
+          vegetation_pars_f%from_file = .FALSE.
+       ENDIF
+
+!
+!--    Read root parameters/distribution and related information
+       IF ( check_existence( var_names, 'soil_pars' ) )  THEN
+          soil_pars_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill,                              &
+                              soil_pars_f%fill,                                &
+                              .FALSE., 'soil_pars' )
+
+          CALL get_attribute( id_surf, char_lod,                               &
+                              soil_pars_f%lod,                                 &
+                              .FALSE., 'soil_pars' )
+
+!
+!--       Inquire number of soil parameters
+          CALL get_dimension_length( id_surf,                                  &
+                                     soil_pars_f%np,                           &
+                                     'nsoil_pars' )
+!
+!--       Read parameters array
+          ALLOCATE( soil_pars_f%pars(0:soil_pars_f%np-1) )
+          CALL get_variable( id_surf, 'nsoil_pars', soil_pars_f%pars )
+
+!
+!--       In case of level of detail 2, also inquire number of vertical
+!--       soil layers, allocate memory and read the respective dimension
+          IF ( soil_pars_f%lod == 2 )  THEN
+             CALL get_dimension_length( id_surf, soil_pars_f%nz, 'zsoil' )
+
+             ALLOCATE( soil_pars_f%layers(0:soil_pars_f%nz-1) )
+             CALL get_variable( id_surf, 'zsoil', soil_pars_f%layers )
+
+          ENDIF
+
+!
+!--       Read soil parameters, depending on level of detail
+          IF ( soil_pars_f%lod == 1 )  THEN
+             ALLOCATE( soil_pars_f%pars_xy(0:soil_pars_f%np-1,                 &
+                                           nys:nyn,nxl:nxr) )
+                  
+             CALL get_variable( id_surf, 'soil_pars', soil_pars_f%pars_xy,     &
+                                nxl, nxr, nys, nyn, 0, soil_pars_f%np-1 )
+
+          ELSEIF ( soil_pars_f%lod == 2 )  THEN
+             ALLOCATE( soil_pars_f%pars_xyz(0:soil_pars_f%np-1,                &
+                                            0:soil_pars_f%nz-1,                &
+                                            nys:nyn,nxl:nxr) )
+             CALL get_variable( id_surf, 'soil_pars',                          &
+                                soil_pars_f%pars_xyz,                          &
+                                nxl, nxr, nys, nyn, 0, soil_pars_f%nz-1,       &
+                                0, soil_pars_f%np-1 )
+
+          ENDIF
+       ELSE
+          soil_pars_f%from_file = .FALSE.
+       ENDIF
+
+!
+!--    Read water parameters and related information
+       IF ( check_existence( var_names, 'water_pars' ) )  THEN
+          water_pars_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill,                              &
+                              water_pars_f%fill,                               &
+                              .FALSE., 'water_pars' )
+!
+!--       Inquire number of water parameters
+          CALL get_dimension_length( id_surf,                                  &
+                                     water_pars_f%np,                          &
+                                     'nwater_pars' )
+!
+!--       Allocate dimension array and input array for water parameters
+          ALLOCATE( water_pars_f%pars(0:water_pars_f%np-1) )
+          ALLOCATE( water_pars_f%pars_xy(0:water_pars_f%np-1,                  &
+                                         nys:nyn,nxl:nxr) )
+!
+!--       Get dimension of water parameters
+          CALL get_variable( id_surf, 'nwater_pars', water_pars_f%pars )
+
+          CALL get_variable( id_surf, 'water_pars', water_pars_f%pars_xy,      &
+                             nxl, nxr, nys, nyn, 0, water_pars_f%np-1 )
+       ELSE
+          water_pars_f%from_file = .FALSE.
+       ENDIF
+!
+!--    Read root area density - parametrized vegetation
+       IF ( check_existence( var_names, 'root_area_dens_s' ) )  THEN
+          root_area_density_lsm_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill,                              &
+                              root_area_density_lsm_f%fill,                    &
+                              .FALSE., 'root_area_dens_s' )
+!
+!--       Obtain number of soil layers from file and allocate variable
+          CALL get_dimension_length( id_surf, root_area_density_lsm_f%nz,      &
+                                     'zsoil' )
+          ALLOCATE( root_area_density_lsm_f%var                                &
+                                        (0:root_area_density_lsm_f%nz-1,       &
+                                         nys:nyn,nxl:nxr) )
+
+!
+!--       Read root-area density
+          CALL get_variable( id_surf, 'root_area_dens_s',                      &
+                             root_area_density_lsm_f%var,                      &
+                             nxl, nxr, nys, nyn,                               &
+                             0, root_area_density_lsm_f%nz-1 )
+
+       ELSE
+          root_area_density_lsm_f%from_file = .FALSE.
+       ENDIF
+!
+!--    Read street type and street crossing
+       IF ( check_existence( var_names, 'street_type' ) )  THEN
+          street_type_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill,                              &
+                              street_type_f%fill, .FALSE.,                     &
+                              'street_type' )
+
+          ALLOCATE ( street_type_f%var(nys:nyn,nxl:nxr)  )
+          
+          CALL get_variable( id_surf, 'street_type', street_type_f%var,        &
+                             nxl, nxr, nys, nyn )
+       ELSE
+          street_type_f%from_file = .FALSE.
+       ENDIF
+
+       IF ( check_existence( var_names, 'street_crossing' ) )  THEN
+          street_crossing_f%from_file = .TRUE.
+          CALL get_attribute( id_surf, char_fill,                              &
+                              street_crossing_f%fill, .FALSE.,                 &
+                              'street_crossing' )
+
+          ALLOCATE ( street_crossing_f%var(nys:nyn,nxl:nxr)  )
+
+          CALL get_variable( id_surf, 'street_crossing',                       &
+                             street_crossing_f%var, nxl, nxr, nys, nyn )
+
+       ELSE
+          street_crossing_f%from_file = .FALSE.
+       ENDIF
+!
+!--    Still missing: root_resolved and building_surface_pars.
+!--    Will be implemented as soon as they are available.
+
+!
+!--    Finally, close input file
+       CALL close_input_file( id_surf )
 #endif
-          ENDIF
-#if defined( __parallel )
-          CALL MPI_BARRIER( comm2d, ierr )
-#endif
-       ENDDO
 !
 !--    End of CPU measurement
@@ -1829,161 +1775,151 @@
        CALL cpu_log( log_point_s(83), 'NetCDF/ASCII input topo', 'start' )
 
-       DO  ii = 0, io_blocks-1
-          IF ( ii == io_group )  THEN
-!
-!--          Input via palm-input data standard
-             IF ( input_pids_static )  THEN
+!
+!--    Input via palm-input data standard
+       IF ( input_pids_static )  THEN
 #if defined ( __netcdf )
 !
-!--             Open file in read-only mode
-                CALL open_read_file( TRIM( input_file_static ) //              &
-                                     TRIM( coupling_char ), id_topo )
-
-!
-!--             At first, inquire all variable names.
-!--             This will be used to check whether an  input variable exist
-!--             or not.
-                CALL inquire_num_variables( id_topo, num_vars )
-!
-!--             Allocate memory to store variable names and inquire them.
-                ALLOCATE( var_names(1:num_vars) )
-                CALL inquire_variable_names( id_topo, var_names )
-!
-!--             Read x, y - dimensions. Only required for consistency checks.
-                CALL get_dimension_length( id_topo, dim_static%nx, 'x' )
-                CALL get_dimension_length( id_topo, dim_static%ny, 'y' )
-                ALLOCATE( dim_static%x(0:dim_static%nx-1) )
-                ALLOCATE( dim_static%y(0:dim_static%ny-1) )
-                CALL get_variable( id_topo, 'x', dim_static%x )
-                CALL get_variable( id_topo, 'y', dim_static%y )
-!
-!--             Terrain height. First, get variable-related _FillValue attribute
-                IF ( check_existence( var_names, 'zt' ) )  THEN
-                   terrain_height_f%from_file = .TRUE.
-                   CALL get_attribute( id_topo, char_fill,                     &
-                                       terrain_height_f%fill,                  &
-                                       .FALSE., 'zt' )
-!
-!--                PE-wise reading of 2D terrain height.
-                   ALLOCATE ( terrain_height_f%var(nys:nyn,nxl:nxr)  )
-                   DO  i = nxl, nxr
-                      CALL get_variable( id_topo, 'zt',                        &
-                                         i, terrain_height_f%var(:,i) )
-                   ENDDO
-                ELSE
-                   terrain_height_f%from_file = .FALSE.
-                ENDIF
-
-!
-!--             Read building height. First, read its _FillValue attribute,
-!--             as well as lod attribute
-                buildings_f%from_file = .FALSE.
-                IF ( check_existence( var_names, 'buildings_2d' ) )  THEN
-                   buildings_f%from_file = .TRUE.
-                   CALL get_attribute( id_topo, char_lod, buildings_f%lod,     &
-                                       .FALSE., 'buildings_2d' )
-
-                   CALL get_attribute( id_topo, char_fill,                     &
-                                       buildings_f%fill1,                      &
-                                       .FALSE., 'buildings_2d' )
-
-!
-!--                Read 2D topography
-                   IF ( buildings_f%lod == 1 )  THEN
-                      ALLOCATE ( buildings_f%var_2d(nys:nyn,nxl:nxr) )
-                      DO  i = nxl, nxr
-                         CALL get_variable( id_topo, 'buildings_2d',           &
-                                            i, buildings_f%var_2d(:,i) )
-                      ENDDO
-                   ELSE
-                      message_string = 'NetCDF attribute lod ' //              &
-                                       '(level of detail) is not set ' //      &
-                                       'properly for buildings_2d.'
-                      CALL message( 'netcdf_data_input_mod', 'PA0457',         &
-                                     1, 2, 0, 6, 0 )
-                   ENDIF
-                ENDIF
-!
-!--             If available, also read 3D building information. If both are
-!--             available, use 3D information.
-                IF ( check_existence( var_names, 'buildings_3d' ) )  THEN
-                   buildings_f%from_file = .TRUE.
-                   CALL get_attribute( id_topo, char_lod, buildings_f%lod,     &
-                                       .FALSE., 'buildings_3d' )
-
-                   CALL get_attribute( id_topo, char_fill,                     &
-                                       buildings_f%fill2,                      &
-                                       .FALSE., 'buildings_3d' )
-
-                   CALL get_dimension_length( id_topo, buildings_f%nz, 'z' )
-
-                   IF ( buildings_f%lod == 2 )  THEN
-                      ALLOCATE( buildings_f%z(nzb:buildings_f%nz-1) )
-                      CALL get_variable( id_topo, 'z', buildings_f%z )
-
-                      ALLOCATE( buildings_f%var_3d(nzb:buildings_f%nz-1,       &
-                                                   nys:nyn,nxl:nxr) )
-                      buildings_f%var_3d = 0
-!
-!--                   Read data PE-wise. Read yz-slices.
-                      DO  i = nxl, nxr
-                         DO  j = nys, nyn
-                            CALL get_variable( id_topo, 'buildings_3d',        &
-                                               i, j,                           &
-                                               buildings_f%var_3d(:,j,i) )
-                         ENDDO
-                      ENDDO
-                   ELSE
-                      message_string = 'NetCDF attribute lod ' //              &
-                                       '(level of detail) is not set ' //      &
-                                       'properly for buildings_3d.'
-                      CALL message( 'netcdf_data_input_mod', 'PA0457',         &
-                                     1, 2, 0, 6, 0 )
-                   ENDIF
-                ENDIF
-!
-!--             Read building IDs and its FillValue attribute. Further required
-!--             for mapping buildings on top of orography.
-                IF ( check_existence( var_names, 'building_id' ) )  THEN
-                   building_id_f%from_file = .TRUE.
-                   CALL get_attribute( id_topo, char_fill,                     &
-                                       building_id_f%fill, .FALSE.,            &
-                                       'building_id' )
-
-
-                   ALLOCATE ( building_id_f%var(nys:nyn,nxl:nxr) )
-                   DO  i = nxl, nxr
-                      CALL get_variable( id_topo, 'building_id',               &
-                                          i, building_id_f%var(:,i) )
-                   ENDDO
-                ELSE
-                   building_id_f%from_file = .FALSE.
-                ENDIF
-!
-!--             Read building_type and required attributes.
-                IF ( check_existence( var_names, 'building_type' ) )  THEN
-                   building_type_f%from_file = .TRUE.
-                   CALL get_attribute( id_topo, char_fill,                     &
-                                       building_type_f%fill, .FALSE.,          &
-                                       'building_type' )
-
-                   ALLOCATE ( building_type_f%var(nys:nyn,nxl:nxr) )
-                   DO  i = nxl, nxr
-                      CALL get_variable( id_topo, 'building_type',             &
-                                         i, building_type_f%var(:,i) )
-                   ENDDO
-                ELSE
-                   building_type_f%from_file = .FALSE.
-                ENDIF
-
-!
-!--             Close topography input file
-                CALL close_input_file( id_topo )
+!--       Open file in read-only mode
+          CALL open_read_file( TRIM( input_file_static ) //                    &
+                               TRIM( coupling_char ), id_topo )
+!
+!--       At first, inquire all variable names.
+!--       This will be used to check whether an  input variable exist
+!--       or not.
+          CALL inquire_num_variables( id_topo, num_vars )
+!
+!--       Allocate memory to store variable names and inquire them.
+          ALLOCATE( var_names(1:num_vars) )
+          CALL inquire_variable_names( id_topo, var_names )
+!
+!--       Read x, y - dimensions. Only required for consistency checks.
+          CALL get_dimension_length( id_topo, dim_static%nx, 'x' )
+          CALL get_dimension_length( id_topo, dim_static%ny, 'y' )
+          ALLOCATE( dim_static%x(0:dim_static%nx-1) )
+          ALLOCATE( dim_static%y(0:dim_static%ny-1) )
+          CALL get_variable( id_topo, 'x', dim_static%x )
+          CALL get_variable( id_topo, 'y', dim_static%y )
+!
+!--       Terrain height. First, get variable-related _FillValue attribute
+          IF ( check_existence( var_names, 'zt' ) )  THEN
+             terrain_height_f%from_file = .TRUE.
+             CALL get_attribute( id_topo, char_fill, terrain_height_f%fill,    &
+                                 .FALSE., 'zt' )
+!
+!--          Input 2D terrain height.
+             ALLOCATE ( terrain_height_f%var(nys:nyn,nxl:nxr)  )
+             
+             CALL get_variable( id_topo, 'zt', terrain_height_f%var,           &
+                                nxl, nxr, nys, nyn )
+
+          ELSE
+             terrain_height_f%from_file = .FALSE.
+          ENDIF
+
+!
+!--       Read building height. First, read its _FillValue attribute,
+!--       as well as lod attribute
+          buildings_f%from_file = .FALSE.
+          IF ( check_existence( var_names, 'buildings_2d' ) )  THEN
+             buildings_f%from_file = .TRUE.
+             CALL get_attribute( id_topo, char_lod, buildings_f%lod,           &
+                                 .FALSE., 'buildings_2d' )
+
+             CALL get_attribute( id_topo, char_fill, buildings_f%fill1,        &
+                                 .FALSE., 'buildings_2d' )
+
+!
+!--          Read 2D buildings
+             IF ( buildings_f%lod == 1 )  THEN
+                ALLOCATE ( buildings_f%var_2d(nys:nyn,nxl:nxr) )
+
+                CALL get_variable( id_topo, 'buildings_2d',                    &
+                                   buildings_f%var_2d,                         &
+                                   nxl, nxr, nys, nyn )
+             ELSE
+                message_string = 'NetCDF attribute lod ' //                    &
+                                 '(level of detail) is not set ' //            &
+                                 'properly for buildings_2d.'
+                CALL message( 'netcdf_data_input_mod', 'NDI000',               &
+                               1, 2, 0, 6, 0 )
+             ENDIF
+          ENDIF
+!
+!--       If available, also read 3D building information. If both are
+!--       available, use 3D information.
+          IF ( check_existence( var_names, 'buildings_3d' ) )  THEN
+             buildings_f%from_file = .TRUE.
+             CALL get_attribute( id_topo, char_lod, buildings_f%lod,           &
+                                 .FALSE., 'buildings_3d' )      
+
+             CALL get_attribute( id_topo, char_fill, buildings_f%fill2,        &
+                                 .FALSE., 'buildings_3d' )
+
+             CALL get_dimension_length( id_topo, buildings_f%nz, 'z' )
+!
+!--          Read 3D buildings
+             IF ( buildings_f%lod == 2 )  THEN
+                ALLOCATE( buildings_f%z(nzb:buildings_f%nz-1) )
+                CALL get_variable( id_topo, 'z', buildings_f%z )
+
+                ALLOCATE( buildings_f%var_3d(nzb:buildings_f%nz-1,             &
+                                             nys:nyn,nxl:nxr) )
+                buildings_f%var_3d = 0
+                
+                CALL get_variable( id_topo, 'buildings_3d',                    &
+                                   buildings_f%var_3d,                         &
+                                   nxl, nxr, nys, nyn, 0, buildings_f%nz-1 )
+             ELSE
+                message_string = 'NetCDF attribute lod ' //                    &
+                                 '(level of detail) is not set ' //            &
+                                 'properly for buildings_3d.'
+                CALL message( 'netcdf_data_input_mod', 'NDI001',               &
+                               1, 2, 0, 6, 0 )
+             ENDIF
+          ENDIF
+!
+!--       Read building IDs and its FillValue attribute. Further required
+!--       for mapping buildings on top of orography.
+          IF ( check_existence( var_names, 'building_id' ) )  THEN
+             building_id_f%from_file = .TRUE.
+             CALL get_attribute( id_topo, char_fill,                           &
+                                 building_id_f%fill, .FALSE.,                  &
+                                 'building_id' )
+
+             ALLOCATE ( building_id_f%var(nys:nyn,nxl:nxr) )
+             
+             CALL get_variable( id_topo, 'building_id', building_id_f%var,     &
+                                nxl, nxr, nys, nyn )
+          ELSE
+             building_id_f%from_file = .FALSE.
+          ENDIF
+!
+!--       Read building_type and required attributes.
+          IF ( check_existence( var_names, 'building_type' ) )  THEN
+             building_type_f%from_file = .TRUE.
+             CALL get_attribute( id_topo, char_fill,                           &
+                                 building_type_f%fill, .FALSE.,                &
+                                 'building_type' )
+
+             ALLOCATE ( building_type_f%var(nys:nyn,nxl:nxr) )
+
+             CALL get_variable( id_topo, 'building_type', building_type_f%var, &
+                                nxl, nxr, nys, nyn )
+
+          ELSE
+             building_type_f%from_file = .FALSE.
+          ENDIF
+!
+!--       Close topography input file
+          CALL close_input_file( id_topo )
 #else
-                CONTINUE
+          CONTINUE
 #endif
 !
-!--          ASCII input
-             ELSEIF ( TRIM( topography ) == 'read_from_file' )  THEN
+!--    ASCII input
+       ELSEIF ( TRIM( topography ) == 'read_from_file' )  THEN
+             
+          DO  ii = 0, io_blocks-1
+             IF ( ii == io_group )  THEN
 
                 OPEN( 90, FILE='TOPOGRAPHY_DATA'//TRIM( coupling_char ),       &
@@ -2021,10 +1957,10 @@
 
              ENDIF
-
-          ENDIF
 #if defined( __parallel )
-          CALL MPI_BARRIER( comm2d, ierr )
+             CALL MPI_BARRIER( comm2d, ierr )
 #endif
-       ENDDO
+          ENDDO
+
+       ENDIF
 !
 !--    End of CPU measurement
@@ -2117,12 +2053,11 @@
        CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names
 
+       LOGICAL      ::  dynamic_3d = .TRUE. !< flag indicating that 3D data is read from dynamic file
+       
        INTEGER(iwp) ::  i          !< running index along x-direction
-       INTEGER(iwp) ::  ii         !< running index for IO blocks
        INTEGER(iwp) ::  id_dynamic !< NetCDF id of dynamic input file
        INTEGER(iwp) ::  j          !< running index along y-direction
        INTEGER(iwp) ::  k          !< running index along z-direction
        INTEGER(iwp) ::  num_vars   !< number of variables in netcdf input file
-       INTEGER(iwp) ::  off_i      !< offset in x-direction used for reading the u-component
-       INTEGER(iwp) ::  off_j      !< offset in y-direction used for reading the v-component
 
        LOGICAL      ::  check_passed !< flag indicating if a check passed
@@ -2148,349 +2083,327 @@
        CALL cpu_log( log_point_s(85), 'NetCDF input init', 'start' )
 
-       DO  ii = 0, io_blocks-1
-          IF ( ii == io_group )  THEN
 #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 of scalar und w grid. Will be used for
-!--          inter- and extrapolation in case of stretched numeric grid.
-!--          This will be removed when Inifor is able to handle stretched grids.
-             CALL get_dimension_length( id_dynamic, init_3d%nzu, 'z'     )
-             CALL get_dimension_length( id_dynamic, init_3d%nzw, 'zw'    )
-             CALL get_dimension_length( id_dynamic, init_3d%nzs, 'zsoil' )
-!
-!--          Read also the horizontal dimensions. These are used just used fo
-!--          checking the compatibility with the PALM grid before reading.
-             CALL get_dimension_length( id_dynamic, init_3d%nx,  'x'  )
-             CALL get_dimension_length( id_dynamic, init_3d%nxu, 'xu' )
-             CALL get_dimension_length( id_dynamic, init_3d%ny,  'y'  )
-             CALL get_dimension_length( id_dynamic, init_3d%nyv, 'yv' )
-
-!
-!--          Check for correct horizontal and vertical dimension. Please note,
-!--          checks are performed directly here and not called from
-!--          check_parameters as some varialbes are still not allocated there.
-!--          Moreover, please note, u- and v-grid has 1 grid point less on
-!--          Inifor grid.
-             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.'
-                CALL message( 'netcdf_data_input_mod', 'NDI003', 1, 2, 0, 6, 0 )
+!--    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 of scalar und w grid. Will be used for
+!--    inter- and extrapolation in case of stretched numeric grid.
+!--    This will be removed when Inifor is able to handle stretched grids.
+       CALL get_dimension_length( id_dynamic, init_3d%nzu, 'z'     )
+       CALL get_dimension_length( id_dynamic, init_3d%nzw, 'zw'    )
+       CALL get_dimension_length( id_dynamic, init_3d%nzs, 'depth' )
+!
+!--    Read also the horizontal dimensions. These are used just used fo
+!--    checking the compatibility with the PALM grid before reading.
+       CALL get_dimension_length( id_dynamic, init_3d%nx,  'x'  )
+       CALL get_dimension_length( id_dynamic, init_3d%nxu, 'xu' )
+       CALL get_dimension_length( id_dynamic, init_3d%ny,  'y'  )
+       CALL get_dimension_length( id_dynamic, init_3d%nyv, 'yv' )
+
+!
+!--    Check for correct horizontal and vertical dimension. Please note,
+!--    checks are performed directly here and not called from
+!--    check_parameters as some varialbes are still not allocated there.
+!--    Moreover, please note, u- and v-grid has 1 grid point less on
+!--    Inifor grid.
+       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.'
+          CALL message( 'netcdf_data_input_mod', 'NDI003', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( init_3d%nzu-1 /= nz )  THEN
+          message_string = 'Number of inifor vertical grid points ' //         &
+                           'does not match the number of numeric grid '//      &
+                           'points.'
+          CALL message( 'netcdf_data_input_mod', 'NDI003', 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, 'z' ) )  THEN
+          ALLOCATE( init_3d%zu_atmos(1:init_3d%nzu) )
+          CALL get_variable( id_dynamic, 'z', init_3d%zu_atmos )
+       ENDIF
+       IF ( check_existence( var_names, 'zw' ) )  THEN
+          ALLOCATE( init_3d%zw_atmos(1:init_3d%nzw) )
+          CALL get_variable( id_dynamic, 'zw', init_3d%zw_atmos )
+       ENDIF
+       IF ( check_existence( var_names, 'depth' ) )  THEN
+          ALLOCATE( init_3d%z_soil(1:init_3d%nzs) )
+          CALL get_variable( id_dynamic, 'depth', init_3d%z_soil )
+       ENDIF
+!
+!--    Read initial geostrophic wind components at t = 0 (index 1 in file).
+!        IF ( check_existence( var_names, 'tend_ug' ) )  THEN
+       IF ( check_existence( var_names, 'ls_forcing_ug' ) )  THEN
+          ALLOCATE( init_3d%ug_init(nzb:nzt+1) )
+!           CALL get_variable_pr( id_dynamic, 'tend_ug', 1,                      &
+!                                 init_3d%ug_init )
+          CALL get_variable_pr( id_dynamic, 'ls_forcing_ug', 1,                &
+                                init_3d%ug_init )
+          init_3d%from_file_ug = .TRUE.
+       ELSE
+          init_3d%from_file_ug = .FALSE.
+       ENDIF
+!        IF ( check_existence( var_names, 'tend_vg' ) )  THEN
+       IF ( check_existence( var_names, 'ls_forcing_vg' ) )  THEN
+          ALLOCATE( init_3d%vg_init(nzb:nzt+1) )
+!           CALL get_variable_pr( id_dynamic, 'tend_vg', 1,                      &
+!                                 init_3d%vg_init )
+          CALL get_variable_pr( id_dynamic, 'ls_forcing_vg', 1,                &
+                                init_3d%vg_init )
+          init_3d%from_file_vg = .TRUE.
+       ELSE
+          init_3d%from_file_vg = .FALSE.
+       ENDIF
+!
+!--    Read inital 3D data of u, v, w, pt and q,
+!--    derived from COSMO model. Read PE-wise yz-slices.
+!--    Please note, the u-, v- and w-component are defined on different
+!--    grids with one element less in the x-, y-,
+!--    and z-direction, respectively. Hence, reading is subdivided
+!--    into separate loops.  
+!--    Read u-component
+       IF ( check_existence( var_names, 'init_u' ) )  THEN
+!
+!--       Read attributes for the fill value and level-of-detail
+          CALL get_attribute( id_dynamic, char_fill, init_3d%fill_u,           &
+                              .FALSE., 'init_u' )
+          CALL get_attribute( id_dynamic, char_lod, init_3d%lod_u,             &
+                              .FALSE., 'init_u' )
+!
+!--       level-of-detail 1 - read initialization profile
+          IF ( init_3d%lod_u == 1 )  THEN
+             ALLOCATE( init_3d%u_init(nzb:nzt+1) )
+             init_3d%u_init = 0.0_wp
+
+             CALL get_variable( id_dynamic, 'init_u',                          &
+                                init_3d%u_init(nzb+1:nzt+1) )
+!
+!--       level-of-detail 2 - read 3D initialization data
+          ELSEIF ( init_3d%lod_u == 2 )  THEN
+
+             CALL get_variable( id_dynamic, 'init_u',                          &
+                                u(nzb+1:nzt+1,nys:nyn,nxlu:nxr),               &
+                                nxlu, nys+1, nzb+1,                            &
+                                nxr-nxlu+1, nyn-nys+1, init_3d%nzu,            &
+                                dynamic_3d )
+          ENDIF
+          init_3d%from_file_u = .TRUE.
+       ENDIF
+!
+!--    Read v-component
+       IF ( check_existence( var_names, 'init_v' ) )  THEN
+!
+!--       Read attributes for the fill value and level-of-detail
+          CALL get_attribute( id_dynamic, char_fill, init_3d%fill_v,           &
+                              .FALSE., 'init_v' )
+          CALL get_attribute( id_dynamic, char_lod, init_3d%lod_v,             &
+                              .FALSE., 'init_v' )
+!
+!--       level-of-detail 1 - read initialization profile
+          IF ( init_3d%lod_v == 1 )  THEN
+             ALLOCATE( init_3d%v_init(nzb:nzt+1) )
+             init_3d%v_init = 0.0_wp
+
+             CALL get_variable( id_dynamic, 'init_v',                          &
+                                init_3d%v_init(nzb+1:nzt+1) )
+
+!
+!--       level-of-detail 2 - read 3D initialization data
+          ELSEIF ( init_3d%lod_v == 2 )  THEN
+
+             CALL get_variable( id_dynamic, 'init_v',                          &
+                                v(nzb+1:nzt+1,nysv:nyn,nxl:nxr),               &
+                                nxl+1, nysv, nzb+1,                            &
+                                nxr-nxl+1, nyn-nysv+1, init_3d%nzu,            &
+                                dynamic_3d )
+             
+          ENDIF
+          init_3d%from_file_v = .TRUE.
+       ENDIF
+!
+!--    Read w-component
+       IF ( check_existence( var_names, 'init_w' ) )  THEN
+!
+!--       Read attributes for the fill value and level-of-detail
+          CALL get_attribute( id_dynamic, char_fill, init_3d%fill_w,           &
+                              .FALSE., 'init_w' )
+          CALL get_attribute( id_dynamic, char_lod, init_3d%lod_w,             &
+                              .FALSE., 'init_w' )
+!
+!--       level-of-detail 1 - read initialization profile
+          IF ( init_3d%lod_w == 1 )  THEN
+             ALLOCATE( init_3d%w_init(nzb:nzt+1) )
+             init_3d%w_init = 0.0_wp
+
+             CALL get_variable( id_dynamic, 'init_w',                          &
+                                init_3d%w_init(nzb+1:nzt) )
+
+!
+!--       level-of-detail 2 - read 3D initialization data
+          ELSEIF ( init_3d%lod_w == 2 )  THEN
+
+             CALL get_variable( id_dynamic, 'init_w',                           &
+                                w(nzb+1:nzt,nys:nyn,nxl:nxr),                   &
+                                nxl+1, nys+1, nzb+1,                            &
+                                nxr-nxl+1, nyn-nys+1, init_3d%nzw,              &
+                                dynamic_3d )
+
+          ENDIF
+          init_3d%from_file_w = .TRUE.
+       ENDIF
+!
+!--    Read potential temperature
+       IF ( .NOT. neutral )  THEN
+          IF ( check_existence( var_names, 'init_pt' ) )  THEN
+!
+!--          Read attributes for the fill value and level-of-detail
+             CALL get_attribute( id_dynamic, char_fill, init_3d%fill_pt,       &
+                                 .FALSE., 'init_pt' )
+             CALL get_attribute( id_dynamic, char_lod, init_3d%lod_pt,         &
+                                 .FALSE., 'init_pt' )
+!
+!--          level-of-detail 1 - read initialization profile
+             IF ( init_3d%lod_pt == 1 )  THEN
+                ALLOCATE( init_3d%pt_init(nzb:nzt+1) )
+
+                CALL get_variable( id_dynamic, 'init_pt',                      &
+                                   init_3d%pt_init(nzb+1:nzt+1) )
+!
+!--             Set Neumann surface boundary condition for initial profil
+                init_3d%pt_init(nzb) = init_3d%pt_init(nzb+1)
+!
+!--          level-of-detail 2 - read 3D initialization data
+             ELSEIF ( init_3d%lod_pt == 2 )  THEN
+
+                CALL get_variable( id_dynamic, 'init_pt',                      &
+                                   pt(nzb+1:nzt+1,nys:nyn,nxl:nxr),            &
+                                   nxl+1, nys+1, nzb+1,                        &
+                                   nxr-nxl+1, nyn-nys+1, init_3d%nzu,          &
+                                   dynamic_3d )
+
+
              ENDIF
-
-             IF ( init_3d%nzu-1 /= nz )  THEN
-                message_string = 'Number of inifor vertical grid points ' //   &
-                                 'does not match the number of numeric grid '//&
-                                 'points.'
-                CALL message( 'netcdf_data_input_mod', 'NDI003', 1, 2, 0, 6, 0 )
+             init_3d%from_file_pt = .TRUE.
+          ENDIF
+       ENDIF
+!
+!--    Read mixing ratio
+       IF ( humidity )  THEN
+          IF ( check_existence( var_names, 'init_qv' ) )  THEN
+!
+!--          Read attributes for the fill value and level-of-detail
+             CALL get_attribute( id_dynamic, char_fill, init_3d%fill_q,        &
+                                 .FALSE., 'init_qv' )
+             CALL get_attribute( id_dynamic, char_lod, init_3d%lod_q,          &
+                                 .FALSE., 'init_qv' )
+!
+!--          level-of-detail 1 - read initialization profile
+             IF ( init_3d%lod_q == 1 )  THEN
+                ALLOCATE( init_3d%q_init(nzb:nzt+1) )
+
+                CALL get_variable( id_dynamic, 'init_qv',                      &
+                                   init_3d%q_init(nzb+1:nzt+1) )
+!
+!--             Set Neumann surface boundary condition for initial profil
+                init_3d%q_init(nzb) = init_3d%q_init(nzb+1)
+
+!
+!--          level-of-detail 2 - read 3D initialization data
+             ELSEIF ( init_3d%lod_q == 2 )  THEN
+             
+                CALL get_variable( id_dynamic, 'init_qv',                      &
+                                   q(nzb+1:nzt+1,nys:nyn,nxl:nxr),             &
+                                   nxl+1, nys+1, nzb+1,                        &
+                                   nxr-nxl+1, nyn-nys+1, init_3d%nzu,          &
+                                   dynamic_3d )
+
+
+
              ENDIF
-!
-!--          Read vertical dimensions. Later, these are required for eventual
-!--          inter- and extrapolations of the initialization data.
-             IF ( check_existence( var_names, 'z' ) )  THEN
-                ALLOCATE( init_3d%zu_atmos(1:init_3d%nzu) )
-                CALL get_variable( id_dynamic, 'z', init_3d%zu_atmos )
+             init_3d%from_file_q = .TRUE.
+          ENDIF
+       ENDIF
+!
+!--    Read soil moisture
+       IF ( land_surface )  THEN
+
+          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_init(0:init_3d%nzs-1) )
+
+                CALL get_variable( id_dynamic, 'init_soil_m',                  &
+                                   init_3d%msoil_init(0:init_3d%nzs-1) )
+!
+!--          level-of-detail 2 - read 3D initialization data
+             ELSEIF ( init_3d%lod_msoil == 2 )  THEN
+                ALLOCATE ( init_3d%msoil(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
+
+               CALL get_variable( id_dynamic, 'init_soil_m',                   &   
+                                  init_3d%msoil(0:init_3d%nzs-1,nys:nyn,nxl:nxr),&
+                                  nxl, nxr, nys, nyn, 0, init_3d%nzs-1 )
+
              ENDIF
-             IF ( check_existence( var_names, 'zw' ) )  THEN
-                ALLOCATE( init_3d%zw_atmos(1:init_3d%nzw) )
-                CALL get_variable( id_dynamic, 'zw', init_3d%zw_atmos )
+             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_init(0:init_3d%nzs-1) )
+
+                CALL get_variable( id_dynamic, 'init_soil_t',                  &
+                                   init_3d%tsoil_init(0:init_3d%nzs-1) )
+
+!
+!--          level-of-detail 2 - read 3D initialization data
+             ELSEIF ( init_3d%lod_tsoil == 2 )  THEN
+                ALLOCATE ( init_3d%tsoil(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
+                
+                CALL get_variable( id_dynamic, 'init_soil_t',                  &   
+                                  init_3d%tsoil(0:init_3d%nzs-1,nys:nyn,nxl:nxr),&
+                                  nxl, nxr, nys, nyn, 0, init_3d%nzs-1 )
              ENDIF
-             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 geostrophic wind components at
-!--          t = 0 (index 1 in file).
-             IF ( check_existence( var_names, 'tend_ug' ) )  THEN
-                ALLOCATE( init_3d%ug_init(nzb:nzt+1) )
-                CALL get_variable_pr( id_dynamic, 'tend_ug', 1,                &
-                                      init_3d%ug_init )
-                init_3d%from_file_ug = .TRUE.
-             ELSE
-                init_3d%from_file_ug = .FALSE.
-             ENDIF
-             IF ( check_existence( var_names, 'tend_vg' ) )  THEN
-                ALLOCATE( init_3d%vg_init(nzb:nzt+1) )
-                CALL get_variable_pr( id_dynamic, 'tend_vg', 1,                &
-                                      init_3d%vg_init )
-                init_3d%from_file_vg = .TRUE.
-             ELSE
-                init_3d%from_file_vg = .FALSE.
-             ENDIF
-!
-!--          Read inital 3D data of u, v, w, pt and q,
-!--          derived from COSMO model. Read PE-wise yz-slices.
-!--          Please note, the u-, v- and w-component are defined on different
-!--          grids with one element less in the x-, y-,
-!--          and z-direction, respectively. Hence, reading is subdivided
-!--          into separate loops. Moreover, i and j are used
-!--          as start index in the NF90 interface.
-!--          The passed arguments for u, and v are (i,j)-1, respectively,
-!--          in contrast to the remaining quantities. This is because in case
-!--          of forcing is applied, the input data for u and v has one
-!--          element less along the x- and y-direction respectively.
-!--          Read u-component
-             IF ( check_existence( var_names, 'init_u' ) )  THEN
-!
-!--             Read attributes for the fill value and level-of-detail
-                CALL get_attribute( id_dynamic, char_fill, init_3d%fill_u,     &
-                                    .FALSE., 'init_u' )
-                CALL get_attribute( id_dynamic, char_lod, init_3d%lod_u,       &
-                                    .FALSE., 'init_u' )
-!
-!--             level-of-detail 1 - read initialization profile
-                IF ( init_3d%lod_u == 1 )  THEN
-                   ALLOCATE( init_3d%u_init(nzb:nzt+1) )
-                   init_3d%u_init = 0.0_wp
-
-                   CALL get_variable( id_dynamic, 'init_u',                    &
-                                      init_3d%u_init(nzb+1:nzt+1) )
-!
-!--             level-of-detail 2 - read 3D initialization data
-                ELSEIF ( init_3d%lod_u == 2 )  THEN
-!
-!--                Set offset value. In case of Dirichlet conditions at the left
-!--                domain boundary, the u component starts at nxl+1. This case,
-!--                the passed start-index for reading the NetCDF data is shifted
-!--                by -1.
-                   off_i = 1 !MERGE( 1, 0, forcing )
-
-                   DO  i = nxlu, nxr
-                      DO  j = nys, nyn
-                         CALL get_variable( id_dynamic, 'init_u', i-off_i, j,  &
-                                            u(nzb+1:nzt+1,j,i) )
-                      ENDDO
-                   ENDDO
-
-                ENDIF
-                init_3d%from_file_u = .TRUE.
-             ENDIF
-!
-!--          Read v-component
-             IF ( check_existence( var_names, 'init_v' ) )  THEN
-!
-!--             Read attributes for the fill value and level-of-detail
-                CALL get_attribute( id_dynamic, char_fill, init_3d%fill_v,     &
-                                    .FALSE., 'init_v' )
-                CALL get_attribute( id_dynamic, char_lod, init_3d%lod_v,       &
-                                    .FALSE., 'init_v' )
-!
-!--             level-of-detail 1 - read initialization profile
-                IF ( init_3d%lod_v == 1 )  THEN
-                   ALLOCATE( init_3d%v_init(nzb:nzt+1) )
-                   init_3d%v_init = 0.0_wp
-
-                   CALL get_variable( id_dynamic, 'init_v',                    &
-                                      init_3d%v_init(nzb+1:nzt+1) )
-
-!
-!--             level-of-detail 2 - read 3D initialization data
-                ELSEIF ( init_3d%lod_v == 2 )  THEN
-!
-!--                Set offset value. In case of Dirichlet conditions at the south
-!--                domain boundary, the v component starts at nys+1. This case,
-!--                the passed start-index for reading the NetCDF data is shifted
-!--                by -1.
-                   off_j = 1 !MERGE( 1, 0, forcing )
-
-                   DO  i = nxl, nxr
-                      DO  j = nysv, nyn
-                         CALL get_variable( id_dynamic, 'init_v', i, j-off_j,  &
-                                            v(nzb+1:nzt+1,j,i) )
-                      ENDDO
-                   ENDDO
-
-                ENDIF
-                init_3d%from_file_v = .TRUE.
-             ENDIF
-!
-!--          Read w-component
-             IF ( check_existence( var_names, 'init_w' ) )  THEN
-!
-!--             Read attributes for the fill value and level-of-detail
-                CALL get_attribute( id_dynamic, char_fill, init_3d%fill_w,     &
-                                    .FALSE., 'init_w' )
-                CALL get_attribute( id_dynamic, char_lod, init_3d%lod_w,       &
-                                    .FALSE., 'init_w' )
-!
-!--             level-of-detail 1 - read initialization profile
-                IF ( init_3d%lod_w == 1 )  THEN
-                   ALLOCATE( init_3d%w_init(nzb:nzt+1) )
-                   init_3d%w_init = 0.0_wp
-
-                   CALL get_variable( id_dynamic, 'init_w',                    &
-                                      init_3d%w_init(nzb+1:nzt) )
-
-!
-!--             level-of-detail 2 - read 3D initialization data
-                ELSEIF ( init_3d%lod_w == 2 )  THEN
-                   DO  i = nxl, nxr
-                      DO  j = nys, nyn
-                         CALL get_variable( id_dynamic, 'init_w', i, j,        &
-                                            w(nzb+1:nzt,j,i) )
-                      ENDDO
-                   ENDDO
-
-                ENDIF
-                init_3d%from_file_w = .TRUE.
-             ENDIF
-!
-!--          Read potential temperature
-             IF ( .NOT. neutral )  THEN
-                IF ( check_existence( var_names, 'init_pt' ) )  THEN
-!
-!--                Read attributes for the fill value and level-of-detail
-                   CALL get_attribute( id_dynamic, char_fill, init_3d%fill_pt, &
-                                       .FALSE., 'init_pt' )
-                   CALL get_attribute( id_dynamic, char_lod, init_3d%lod_pt,   &
-                                       .FALSE., 'init_pt' )
-!
-!--                level-of-detail 1 - read initialization profile
-                   IF ( init_3d%lod_pt == 1 )  THEN
-                      ALLOCATE( init_3d%pt_init(nzb:nzt+1) )
-
-                      CALL get_variable( id_dynamic, 'init_pt',                &
-                                         init_3d%pt_init(nzb+1:nzt+1) )
-!
-!--                   Set Neumann surface boundary condition for initial profil
-                      init_3d%pt_init(nzb) = init_3d%pt_init(nzb+1)
-!
-!--                level-of-detail 2 - read 3D initialization data
-                   ELSEIF ( init_3d%lod_pt == 2 )  THEN
-                      DO  i = nxl, nxr
-                         DO  j = nys, nyn
-                            CALL get_variable( id_dynamic, 'init_pt', i, j,    &
-                                               pt(nzb+1:nzt+1,j,i) )
-                         ENDDO
-                      ENDDO
-
-                   ENDIF
-                   init_3d%from_file_pt = .TRUE.
-                ENDIF
-             ENDIF
-!
-!--          Read mixing ratio
-             IF ( humidity )  THEN
-                IF ( check_existence( var_names, 'init_qv' ) )  THEN
-!
-!--                Read attributes for the fill value and level-of-detail
-                   CALL get_attribute( id_dynamic, char_fill, init_3d%fill_q,  &
-                                       .FALSE., 'init_qv' )
-                   CALL get_attribute( id_dynamic, char_lod, init_3d%lod_q,    &
-                                       .FALSE., 'init_qv' )
-!
-!--                level-of-detail 1 - read initialization profile
-                   IF ( init_3d%lod_q == 1 )  THEN
-                      ALLOCATE( init_3d%q_init(nzb:nzt+1) )
-
-                      CALL get_variable( id_dynamic, 'init_qv',               &
-                                         init_3d%q_init(nzb+1:nzt+1) )
-!
-!--                   Set Neumann surface boundary condition for initial profil
-                      init_3d%q_init(nzb) = init_3d%q_init(nzb+1)
-
-!
-!--                level-of-detail 2 - read 3D initialization data
-                   ELSEIF ( init_3d%lod_q == 2 )  THEN
-                      DO  i = nxl, nxr
-                         DO  j = nys, nyn
-                            CALL get_variable( id_dynamic, 'init_qv', i, j,    &
-                                               q(nzb+1:nzt+1,j,i) )
-                         ENDDO
-                      ENDDO
-
-                   ENDIF
-                   init_3d%from_file_q = .TRUE.
-                ENDIF
-             ENDIF
-!
-!--          Read soil moisture
-             IF ( land_surface )  THEN
-
-                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_init(0:init_3d%nzs-1) )
-
-                      CALL get_variable( id_dynamic, 'init_soil_m',            &
-                                         init_3d%msoil_init(0:init_3d%nzs-1) )
-!
-!--                level-of-detail 2 - read 3D initialization data
-                   ELSEIF ( init_3d%lod_msoil == 2 )  THEN
-                      ALLOCATE ( init_3d%msoil(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
-                      DO  i = nxl, nxr
-                         DO  j = nys, nyn
-                            CALL get_variable( id_dynamic, 'init_soil_m', i, j,&
-                                               init_3d%msoil(0:init_3d%nzs-1,j,i) )
-                         ENDDO
-                      ENDDO
-                   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_init(0:init_3d%nzs-1) )
-
-                      CALL get_variable( id_dynamic, 'init_soil_t',            &
-                                         init_3d%tsoil_init(0:init_3d%nzs-1) )
-
-!
-!--                level-of-detail 2 - read 3D initialization data
-                   ELSEIF ( init_3d%lod_tsoil == 2 )  THEN
-                      ALLOCATE ( init_3d%tsoil(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
-                      DO  i = nxl, nxr
-                         DO  j = nys, nyn
-                            CALL get_variable( id_dynamic, 'init_soil_t', i, j,&
-                                               init_3d%tsoil(0:init_3d%nzs-1,j,i) )
-                         ENDDO
-                      ENDDO
-                   ENDIF
-                   init_3d%from_file_tsoil = .TRUE.
-                ENDIF
-             ENDIF
-!
-!--          Close input file
-             CALL close_input_file( id_dynamic )
+             init_3d%from_file_tsoil = .TRUE.
+          ENDIF
+       ENDIF
+!
+!--    Close input file
+       CALL close_input_file( id_dynamic )
 #endif
-          ENDIF
-#if defined( __parallel )
-          CALL MPI_BARRIER( comm2d, ierr )
-#endif
-       ENDDO
 !
 !--    End of CPU measurement
@@ -2604,7 +2517,7 @@
        IMPLICIT NONE
 
-
+       LOGICAL      ::  dynamic_3d = .TRUE. !< flag indicating that 3D data is read from dynamic file
+       
        INTEGER(iwp) ::  i          !< running index along x-direction
-       INTEGER(iwp) ::  ii         !< running index for IO blocks
        INTEGER(iwp) ::  id_dynamic !< NetCDF id of dynamic input file
        INTEGER(iwp) ::  j          !< running index along y-direction
@@ -2624,368 +2537,238 @@
        CALL cpu_log( log_point_s(86), 'NetCDF input forcing', 'start' )
 
-       DO  ii = 0, io_blocks-1
-          IF ( ii == io_group )  THEN
 #if defined ( __netcdf )
 !
-!--          Open file in read-only mode
-             CALL open_read_file( TRIM( input_file_dynamic ) //                &
-                                  TRIM( coupling_char ), id_dynamic )
-!
-!--          Initialize INIFOR forcing.
-             IF ( .NOT. force%init )  THEN
-!
-!--             At first, inquire all variable names.
-                CALL inquire_num_variables( id_dynamic, num_vars )
-!
-!--             Allocate memory to store variable names.
-                ALLOCATE( force%var_names(1:num_vars) )
-                CALL inquire_variable_names( id_dynamic, force%var_names )
-!
-!--             Read time dimension, allocate memory and finally read time array
-                CALL get_dimension_length( id_dynamic, force%nt, 'time' )
-
-                IF ( check_existence( force%var_names, 'time' ) )  THEN
-                   ALLOCATE( force%time(0:force%nt-1) )
-                   CALL get_variable( id_dynamic, 'time', force%time )
-                ENDIF
-!
-!--             Read vertical dimension of scalar und w grid
-                CALL get_dimension_length( id_dynamic, force%nzu, 'z' )
-                CALL get_dimension_length( id_dynamic, force%nzw, 'zw' )
-
-                IF ( check_existence( force%var_names, 'z' ) )  THEN
-                   ALLOCATE( force%zu_atmos(1:force%nzu) )
-                   CALL get_variable( id_dynamic, 'z', force%zu_atmos )
-                ENDIF
-                IF ( check_existence( force%var_names, 'zw' ) )  THEN
-                   ALLOCATE( force%zw_atmos(1:force%nzw) )
-                   CALL get_variable( id_dynamic, 'zw', force%zw_atmos )
-                ENDIF
-
-!
-!--             Read surface pressure
-                IF ( check_existence( force%var_names,                         &
-                                  'surface_forcing_surface_pressure' ) )  THEN
-                   ALLOCATE( force%surface_pressure(0:force%nt-1) )
-                   CALL get_variable( id_dynamic,                              &
-                                      'surface_forcing_surface_pressure',      &
-                                      force%surface_pressure )
-                ENDIF
-!
-!--             Set control flag to indicate that initialization is already done
-                force%init = .TRUE.
-
-             ENDIF
-
-!
-!--          Obtain time index for current input starting at 0.
-!--          @todo: At the moment time, in INIFOR and simulated time correspond
-!--                 to each other. If required, adjust to daytime.
-             force%tind = MINLOC( ABS( force%time - simulated_time ), DIM = 1 )&
-                          - 1
-             force%tind_p = force%tind + 1
-!
-!--          Read geostrophic wind components. In case of forcing, this is only
-!--          required if cyclic boundary conditions are applied.
-             IF ( bc_lr_cyc  .AND.  bc_ns_cyc )  THEN
-                DO  t = force%tind, force%tind_p
-                   CALL get_variable_pr( id_dynamic, 'tend_ug', t+1,           &
-                                         force%ug(t-force%tind,:) )
-                   CALL get_variable_pr( id_dynamic, 'tend_vg', t+1,           &
-                                         force%ug(t-force%tind,:) )
-                ENDDO
-             ENDIF
-!
-!--          Read data at lateral and top boundaries. Please note, at left and
-!--          right domain boundary, yz-layers are read for u, v, w, pt and q.
-!--          For the v-component, the data starts at nysv, while for the other
-!--          quantities the data starts at nys. This is equivalent at the north
-!--          and south domain boundary for the u-component.
-!--          The function get_variable_bc assumes the start indices with respect
-!--          to the netcdf file convention (data starts at index 1). For this
-!--          reason, nys+1 / nxl+1 are passed instead of nys / nxl. For the
-!--          the u- and v-component at the north/south, and left/right boundary,
-!--          nxlu and nysv are passed, respectively, since these always starts
-!--          at index 1 in case of forcing.
-
-             IF ( force_bound_l )  THEN
-                DO  j = nys, nyn
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic, 'ls_forcing_left_u',   &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      j+1, 1,                                  &
-                                      force%u_left(t-force%tind,nzb+1:nzt+1,j) )
-                   ENDDO
-                ENDDO
-
-                DO  j = nysv, nyn
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic, 'ls_forcing_left_v',   &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      j, 1,                                    &
-                                      force%v_left(t-force%tind,nzb+1:nzt+1,j) )
-                   ENDDO
-                ENDDO
-                DO  j = nys, nyn
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic,                        &
-                                      'ls_forcing_left_w',                     &
-                                      t+1,                                     &
-                                      nzb+1, nzt-(nzb+1) + 1,                  &
-                                      j+1, 1,                                  &
-                                      force%w_left(t-force%tind,nzb+1:nzt,j) )
-                   ENDDO
-                ENDDO
-                IF ( .NOT. neutral )  THEN
-                   DO  j = nys, nyn
-                      DO  t = force%tind, force%tind_p
-                         CALL get_variable_bc( id_dynamic,                     &
-                                      'ls_forcing_left_pt',                    &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      j+1, 1,                                  &
-                                      force%pt_left(t-force%tind,nzb+1:nzt+1,j) )
-                      ENDDO
-                   ENDDO
-                ENDIF
-                IF ( humidity )  THEN
-                   DO  j = nys, nyn
-                      DO  t = force%tind, force%tind_p
-                         CALL get_variable_bc( id_dynamic,                     &
-                                      'ls_forcing_left_qv',                    &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      j+1, 1,                                  &
-                                      force%q_left(t-force%tind,nzb+1:nzt+1,j) )
-                      ENDDO
-                   ENDDO
-                ENDIF
-             ENDIF
-
-             IF ( force_bound_r )  THEN
-                DO  j = nys, nyn
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic, 'ls_forcing_right_u',  &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      j+1, 1,                                  &
-                                      force%u_right(t-force%tind,nzb+1:nzt+1,j) )
-                   ENDDO
-                ENDDO
-                DO  j = nysv, nyn
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic, 'ls_forcing_right_v',  &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      j, 1,                                    &
-                                      force%v_right(t-force%tind,nzb+1:nzt+1,j) )
-                   ENDDO
-                ENDDO
-                DO  j = nys, nyn
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic, 'ls_forcing_right_w',  &
-                                      t+1,                                     &
-                                      nzb+1, nzt-(nzb+1)+1,                    &
-                                      j+1, 1,                                  &
-                                      force%w_right(t-force%tind,nzb+1:nzt,j) )
-                   ENDDO
-                ENDDO
-                IF ( .NOT. neutral )  THEN
-                   DO  j = nys, nyn
-                      DO  t = force%tind, force%tind_p
-                         CALL get_variable_bc( id_dynamic,                     &
-                                      'ls_forcing_right_pt',                   &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      j+1, 1,                                  &
-                                      force%pt_right(t-force%tind,nzb+1:nzt+1,j) )
-                      ENDDO
-                   ENDDO
-                ENDIF
-                IF ( humidity )  THEN
-                   DO  j = nys, nyn
-                      DO  t = force%tind, force%tind_p
-                         CALL get_variable_bc( id_dynamic,                     &
-                                      'ls_forcing_right_qv',                   &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      j+1, 1,                                  &
-                                      force%q_right(t-force%tind,nzb+1:nzt+1,j) )
-                      ENDDO
-                   ENDDO
-                ENDIF
-             ENDIF
-
-             IF ( force_bound_n )  THEN
-                DO  i = nxlu, nxr
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic, 'ls_forcing_north_u',  &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      i, 1,                                    &
-                                      force%u_north(t-force%tind,nzb+1:nzt+1,i) )
-                   ENDDO
-                ENDDO
-                DO  i = nxl, nxr
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic, 'ls_forcing_north_v',  &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      i+1, 1,                                  &
-                                      force%v_north(t-force%tind,nzb+1:nzt+1,i) )
-                   ENDDO
-                ENDDO
-                DO  i = nxl, nxr
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic, 'ls_forcing_north_w',  &
-                                      t+1,                                     &
-                                      nzb+1, nzt-(nzb+1)+1,                    &
-                                      i+1, 1,                                  &
-                                      force%w_north(t-force%tind,nzb+1:nzt,i) )
-                   ENDDO
-                ENDDO
-                IF ( .NOT. neutral )  THEN
-                   DO  i = nxl, nxr
-                      DO  t = force%tind, force%tind_p
-                         CALL get_variable_bc( id_dynamic,                     &
-                                      'ls_forcing_north_pt',                   &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      i+1, 1,                                  &
-                                      force%pt_north(t-force%tind,nzb+1:nzt+1,i) )
-                      ENDDO
-                   ENDDO
-                ENDIF
-                IF ( humidity )  THEN
-                   DO  i = nxl, nxr
-                      DO  t = force%tind, force%tind_p
-                         CALL get_variable_bc( id_dynamic,                     &
-                                      'ls_forcing_north_qv',                   &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      i+1, 1,                                  &
-                                      force%q_north(t-force%tind,nzb+1:nzt+1,i) )
-                      ENDDO
-                   ENDDO
-                ENDIF
-             ENDIF
-
-             IF ( force_bound_s )  THEN
-                DO  i = nxlu, nxr
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic, 'ls_forcing_south_u',  &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      i, 1,                                    &
-                                      force%u_south(t-force%tind,nzb+1:nzt+1,i) )
-                   ENDDO
-                ENDDO
-                DO  i = nxl, nxr
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic, 'ls_forcing_south_v',  &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      i+1, 1,                                  &
-                                      force%v_south(t-force%tind,nzb+1:nzt+1,i) )
-                   ENDDO
-                ENDDO
-                DO  i = nxl, nxr
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic, 'ls_forcing_south_w',  &
-                                      t+1,                                     &
-                                      nzb+1, nzt-(nzb+1)+1,                    &
-                                      i+1, 1,                                  &
-                                      force%w_south(t-force%tind,nzb+1:nzt,i) )
-                   ENDDO
-                ENDDO
-                IF ( .NOT. neutral )  THEN
-                   DO  i = nxl, nxr
-                      DO  t = force%tind, force%tind_p
-                         CALL get_variable_bc( id_dynamic,                     &
-                                      'ls_forcing_south_pt',                   &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      i+1, 1,                                  &
-                                      force%pt_south(t-force%tind,nzb+1:nzt+1,i) )
-                      ENDDO
-                   ENDDO
-                ENDIF
-                IF ( humidity )  THEN
-                   DO  i = nxl, nxr
-                      DO  t = force%tind, force%tind_p
-                         CALL get_variable_bc( id_dynamic,                     &
-                                      'ls_forcing_south_qv',                   &
-                                      t+1,                                     &
-                                      nzb+1, nzt+1-(nzb+1)+1,                  &
-                                      i+1, 1,                                  &
-                                      force%q_south(t-force%tind,nzb+1:nzt+1,i) )
-                      ENDDO
-                   ENDDO
-                ENDIF
-             ENDIF
-!
-!--          Top boundary
-             DO  i = nxlu, nxr
-                DO  t = force%tind, force%tind_p
-                   CALL get_variable_bc( id_dynamic, 'ls_forcing_top_u',       &
-                                   t+1,                                        &
-                                   nys+1, nyn-nys+1,                           &
-                                   i, 1,                                       &
-                                   force%u_top(t-force%tind,nys:nyn,i) )
-                ENDDO
-             ENDDO
-             DO  i = nxl, nxr
-                DO  t = force%tind, force%tind_p
-                   CALL get_variable_bc( id_dynamic, 'ls_forcing_top_v',       &
-                                   t+1,                                        &
-                                   nysv, nyn-nysv+1,                           &
-                                   i+1, 1,                                     &
-                                   force%v_top(t-force%tind,nysv:nyn,i) )
-                ENDDO
-             ENDDO
-             DO  i = nxl, nxr
-                DO  t = force%tind, force%tind_p
-                   CALL get_variable_bc( id_dynamic, 'ls_forcing_top_w',       &
-                                   t+1,                                        &
-                                   nys+1, nyn-nys+1,                           &
-                                   i+1, 1,                                     &
-                                   force%w_top(t-force%tind,nys:nyn,i) )
-                ENDDO
-             ENDDO
-             IF ( .NOT. neutral )  THEN
-                DO  i = nxl, nxr
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic, 'ls_forcing_top_pt',   &
-                                      t+1,                                     &
-                                      nys+1, nyn-nys+1,                        &
-                                      i+1, 1,                                  &
-                                      force%pt_top(t-force%tind,nys:nyn,i) )
-                   ENDDO
-                ENDDO
-             ENDIF
-             IF ( humidity )  THEN
-                DO  i = nxl, nxr
-                   DO  t = force%tind, force%tind_p
-                      CALL get_variable_bc( id_dynamic, 'ls_forcing_top_qv',   &
-                                      t+1,                                     &
-                                      nys+1, nyn-nys+1,                        &
-                                      i+1, 1,                                  &
-                                      force%q_top(t-force%tind,nys:nyn,i) )
-                   ENDDO
-                ENDDO
-             ENDIF
-
-!
-!--          Close input file
-             CALL close_input_file( id_dynamic )
+!--    Open file in read-only mode
+       CALL open_read_file( TRIM( input_file_dynamic ) //                      &
+                            TRIM( coupling_char ), id_dynamic )
+!
+!--    Initialize INIFOR forcing.
+       IF ( .NOT. force%init )  THEN
+!
+!--       At first, inquire all variable names.
+          CALL inquire_num_variables( id_dynamic, num_vars )
+!
+!--       Allocate memory to store variable names.
+          ALLOCATE( force%var_names(1:num_vars) )
+          CALL inquire_variable_names( id_dynamic, force%var_names )
+!
+!--       Read time dimension, allocate memory and finally read time array
+          CALL get_dimension_length( id_dynamic, force%nt, 'time' )
+
+          IF ( check_existence( force%var_names, 'time' ) )  THEN
+             ALLOCATE( force%time(0:force%nt-1) )
+             CALL get_variable( id_dynamic, 'time', force%time )
+          ENDIF
+!
+!--       Read vertical dimension of scalar und w grid
+          CALL get_dimension_length( id_dynamic, force%nzu, 'z' )
+          CALL get_dimension_length( id_dynamic, force%nzw, 'zw' )
+
+          IF ( check_existence( force%var_names, 'z' ) )  THEN
+             ALLOCATE( force%zu_atmos(1:force%nzu) )
+             CALL get_variable( id_dynamic, 'z', force%zu_atmos )
+          ENDIF
+          IF ( check_existence( force%var_names, 'zw' ) )  THEN
+             ALLOCATE( force%zw_atmos(1:force%nzw) )
+             CALL get_variable( id_dynamic, 'zw', force%zw_atmos )
+          ENDIF
+
+!
+!--       Read surface pressure
+          IF ( check_existence( force%var_names,                               &
+                            'surface_forcing_surface_pressure' ) )  THEN
+             ALLOCATE( force%surface_pressure(0:force%nt-1) )
+             CALL get_variable( id_dynamic,                                    &
+                                'surface_forcing_surface_pressure',            &
+                                force%surface_pressure )
+          ENDIF
+!
+!--       Set control flag to indicate that initialization is already done
+          force%init = .TRUE.
+
+       ENDIF
+
+!
+!--    Obtain time index for current input starting at 0.
+!--    @todo: At the moment time, in INIFOR and simulated time correspond
+!--           to each other. If required, adjust to daytime.
+       force%tind = MINLOC( ABS( force%time - simulated_time ), DIM = 1 )      &
+                    - 1
+       force%tind_p = force%tind + 1       
+!
+!--    Read geostrophic wind components. In case of forcing, this is only
+!--    required if cyclic boundary conditions are applied.
+       IF ( bc_lr_cyc  .AND.  bc_ns_cyc )  THEN
+          DO  t = force%tind, force%tind_p
+!              CALL get_variable_pr( id_dynamic, 'tend_ug', t+1,           &
+!                                    force%ug(t-force%tind,:) )
+!              CALL get_variable_pr( id_dynamic, 'tend_vg', t+1,           &
+!                                    force%ug(t-force%tind,:) )
+             CALL get_variable_pr( id_dynamic, 'ls_forcing_ug', t+1,           &
+                                   force%ug(t-force%tind,:) )
+             CALL get_variable_pr( id_dynamic, 'ls_forcing_vg', t+1,           &
+                                   force%ug(t-force%tind,:) )
+          ENDDO
+       ENDIF
+!
+!--    Read data at lateral and top boundaries. Please note, at left and
+!--    right domain boundary, yz-layers are read for u, v, w, pt and q.
+!--    For the v-component, the data starts at nysv, while for the other
+!--    quantities the data starts at nys. This is equivalent at the north
+!--    and south domain boundary for the u-component.
+       IF ( force_bound_l )  THEN
+          CALL get_variable( id_dynamic, 'ls_forcing_left_u',                  &
+                           force%u_left(0:1,nzb+1:nzt+1,nys:nyn),              &
+                           nys+1, nzb+1, force%tind+1,                         &
+                           nyn-nys+1, force%nzu, 2, dynamic_3d )
+          
+          CALL get_variable( id_dynamic, 'ls_forcing_left_v',                  &
+                           force%v_left(0:1,nzb+1:nzt+1,nysv:nyn),             &
+                           nysv, nzb+1, force%tind+1,                          &
+                           nyn-nysv+1, force%nzu, 2, dynamic_3d )
+
+          CALL get_variable( id_dynamic, 'ls_forcing_left_w',                  &
+                           force%w_left(0:1,nzb+1:nzt,nys:nyn),                &
+                           nys+1, nzb+1, force%tind+1,                         &
+                           nyn-nys+1, force%nzw, 2, dynamic_3d )
+
+          IF ( .NOT. neutral )  THEN
+             CALL get_variable( id_dynamic, 'ls_forcing_left_pt',              &
+                           force%pt_left(0:1,nzb+1:nzt+1,nys:nyn),             &
+                           nys+1, nzb+1, force%tind+1,                         &
+                           nyn-nys+1, force%nzu, 2, dynamic_3d )
+          ENDIF
+          IF ( humidity )  THEN
+             CALL get_variable( id_dynamic, 'ls_forcing_left_qv',              &
+                           force%q_left(0:1,nzb+1:nzt+1,nys:nyn),              &
+                           nys+1, nzb+1, force%tind+1,                         &
+                           nyn-nys+1, force%nzu, 2, dynamic_3d )
+          ENDIF
+       ENDIF
+
+       IF ( force_bound_r )  THEN
+          CALL get_variable( id_dynamic, 'ls_forcing_right_u',                 &
+                           force%u_right(0:1,nzb+1:nzt+1,nys:nyn),             &
+                           nys+1, nzb+1, force%tind+1,                         &
+                           nyn-nys+1, force%nzu, 2, dynamic_3d )
+                           
+          CALL get_variable( id_dynamic, 'ls_forcing_right_v',                 &
+                           force%v_right(0:1,nzb+1:nzt+1,nysv:nyn),            &
+                           nysv, nzb+1, force%tind+1,                          &
+                           nyn-nysv+1, force%nzu, 2, dynamic_3d )
+                           
+          CALL get_variable( id_dynamic, 'ls_forcing_right_w',                 &
+                           force%w_right(0:1,nzb+1:nzt,nys:nyn),               &
+                           nys+1, nzb+1, force%tind+1,                         &
+                           nyn-nys+1, force%nzw, 2, dynamic_3d )
+                           
+          IF ( .NOT. neutral )  THEN
+             CALL get_variable( id_dynamic, 'ls_forcing_right_pt',             &
+                           force%pt_right(0:1,nzb+1:nzt+1,nys:nyn),            &
+                           nys+1, nzb+1, force%tind+1,                         &
+                           nyn-nys+1, force%nzu, 2, dynamic_3d )
+          ENDIF
+          IF ( humidity )  THEN
+             CALL get_variable( id_dynamic, 'ls_forcing_right_qv',             &
+                           force%q_right(0:1,nzb+1:nzt+1,nys:nyn),             &
+                           nys+1, nzb+1, force%tind+1,                         &
+                           nyn-nys+1, force%nzu, 2, dynamic_3d )
+          ENDIF
+       ENDIF
+
+       IF ( force_bound_n )  THEN
+       
+          CALL get_variable( id_dynamic, 'ls_forcing_north_u',                 &
+                           force%u_north(0:1,nzb+1:nzt+1,nxlu:nxr),            &
+                           nxlu, nzb+1, force%tind+1,                          &
+                           nxr-nxlu+1, force%nzu, 2, dynamic_3d )
+
+          CALL get_variable( id_dynamic, 'ls_forcing_north_v',                 &
+                           force%v_north(0:1,nzb+1:nzt+1,nxl:nxr),             &
+                           nxl+1, nzb+1, force%tind+1,                         &
+                           nxr-nxl+1, force%nzu, 2, dynamic_3d )
+                           
+          CALL get_variable( id_dynamic, 'ls_forcing_north_w',                 &
+                           force%w_north(0:1,nzb+1:nzt,nxl:nxr),               &
+                           nxl+1, nzb+1, force%tind+1,                         &
+                           nxr-nxl+1, force%nzw, 2, dynamic_3d )
+                           
+          IF ( .NOT. neutral )  THEN
+             CALL get_variable( id_dynamic, 'ls_forcing_north_pt',             &
+                           force%pt_north(0:1,nzb+1:nzt+1,nxl:nxr),            &
+                           nxl+1, nzb+1, force%tind+1,                         &
+                           nxr-nxl+1, force%nzu, 2, dynamic_3d )
+          ENDIF
+          IF ( humidity )  THEN
+             CALL get_variable( id_dynamic, 'ls_forcing_north_qv',             &
+                           force%q_north(0:1,nzb+1:nzt+1,nxl:nxr),             &
+                           nxl+1, nzb+1, force%tind+1,                         &
+                           nxr-nxl+1, force%nzu, 2, dynamic_3d )
+          ENDIF
+       ENDIF
+
+       IF ( force_bound_s )  THEN
+          CALL get_variable( id_dynamic, 'ls_forcing_south_u',                 &
+                           force%u_south(0:1,nzb+1:nzt+1,nxlu:nxr),            &
+                           nxlu, nzb+1, force%tind+1,                          &
+                           nxr-nxlu+1, force%nzu, 2, dynamic_3d )
+
+          CALL get_variable( id_dynamic, 'ls_forcing_south_v',                 &
+                           force%v_south(0:1,nzb+1:nzt+1,nxl:nxr),             &
+                           nxl+1, nzb+1, force%tind+1,                         &
+                           nxr-nxl+1, force%nzu, 2, dynamic_3d )
+                           
+          CALL get_variable( id_dynamic, 'ls_forcing_south_w',                 &
+                           force%w_south(0:1,nzb+1:nzt,nxl:nxr),               &
+                           nxl+1, nzb+1, force%tind+1,                         &
+                           nxr-nxl+1, force%nzw, 2, dynamic_3d )
+                           
+          IF ( .NOT. neutral )  THEN
+             CALL get_variable( id_dynamic, 'ls_forcing_south_pt',             &
+                           force%pt_south(0:1,nzb+1:nzt+1,nxl:nxr),            &
+                           nxl+1, nzb+1, force%tind+1,                         &
+                           nxr-nxl+1, force%nzu, 2, dynamic_3d )
+          ENDIF
+          IF ( humidity )  THEN
+             CALL get_variable( id_dynamic, 'ls_forcing_south_qv',             &
+                           force%q_south(0:1,nzb+1:nzt+1,nxl:nxr),             &
+                           nxl+1, nzb+1, force%tind+1,                         &
+                           nxr-nxl+1, force%nzu, 2, dynamic_3d )
+          ENDIF
+       ENDIF
+!
+!--    Top boundary
+       CALL get_variable( id_dynamic, 'ls_forcing_top_u',                      &
+                             force%u_top(0:1,nys:nyn,nxlu:nxr),                &
+                             nxlu, nys+1, force%tind+1,                        &
+                             nxr-nxlu+1, nyn-nys+1, 2, dynamic_3d )
+
+       CALL get_variable( id_dynamic, 'ls_forcing_top_v',                      &
+                             force%v_top(0:1,nysv:nyn,nxl:nxr),                &
+                             nxl+1, nysv, force%tind+1,                        &
+                             nxr-nxl+1, nyn-nysv+1, 2, dynamic_3d )
+                             
+       CALL get_variable( id_dynamic, 'ls_forcing_top_w',                      &
+                             force%w_top(0:1,nys:nyn,nxl:nxr),                 &
+                             nxl+1, nys+1, force%tind+1,                       &
+                             nxr-nxl+1, nyn-nys+1, 2, dynamic_3d )
+                             
+       IF ( .NOT. neutral )  THEN
+          CALL get_variable( id_dynamic, 'ls_forcing_top_pt',                  &
+                                force%pt_top(0:1,nys:nyn,nxl:nxr),             &
+                                nxl+1, nys+1, force%tind+1,                    &
+                                nxr-nxl+1, nyn-nys+1, 2, dynamic_3d )
+       ENDIF
+       IF ( humidity )  THEN
+          CALL get_variable( id_dynamic, 'ls_forcing_top_qv',                  &
+                                force%q_top(0:1,nys:nyn,nxl:nxr),              &
+                                nxl+1, nys+1, force%tind+1,                    &
+                                nxr-nxl+1, nyn-nys+1, 2, dynamic_3d )
+       ENDIF
+
+!
+!--    Close input file
+       CALL close_input_file( id_dynamic )
 #endif
-          ENDIF
-#if defined( __parallel )
-          CALL MPI_BARRIER( comm2d, ierr )
-#endif
-       ENDDO
-
 !
 !--    End of CPU measurement
@@ -3874,5 +3657,5 @@
 
        nc_stat = NF90_CLOSE( id )
-       CALL handle_error( 'close', 537 )
+       CALL handle_error( 'close', 540 )
 #endif
     END SUBROUTINE close_input_file
@@ -3901,8 +3684,6 @@
            nc_stat = NF90_OPEN( filename, NF90_NOWRITE, id )
            collective_read = .FALSE.
-           write(9,*) 'open ',TRIM(filename),' as serial NetCDF file on every PE'
        ELSE
            collective_read = .TRUE.
-           write(9,*) 'open ',TRIM(filename),' using parallel NetCDF'
        END IF
 #else
@@ -3911,5 +3692,5 @@
 #endif
 
-       CALL handle_error( 'open_read_file', 536 )
+       CALL handle_error( 'open_read_file', 539 )
 
 #endif
@@ -3942,5 +3723,5 @@
        IF ( global )  THEN
           nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
-          CALL handle_error( 'get_attribute_int32 global', 522 )
+          CALL handle_error( 'get_attribute_int32 global', 522, attribute_name )
 !
 !--    Read attributes referring to a single variable. Therefore, first inquire
@@ -3948,7 +3729,7 @@
        ELSE
           nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
-          CALL handle_error( 'get_attribute_int32', 522 )
+          CALL handle_error( 'get_attribute_int32', 522, attribute_name )
           nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
-          CALL handle_error( 'get_attribute_int32', 522 )
+          CALL handle_error( 'get_attribute_int32', 522, attribute_name )
        ENDIF
 #endif
@@ -3981,5 +3762,5 @@
        IF ( global )  THEN
           nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
-          CALL handle_error( 'get_attribute_int8 global', 523 )
+          CALL handle_error( 'get_attribute_int8 global', 523, attribute_name )
 !
 !--    Read attributes referring to a single variable. Therefore, first inquire
@@ -3987,7 +3768,7 @@
        ELSE
           nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
-          CALL handle_error( 'get_attribute_int8', 523 )
+          CALL handle_error( 'get_attribute_int8', 523, attribute_name )
           nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
-          CALL handle_error( 'get_attribute_int8', 523 )
+          CALL handle_error( 'get_attribute_int8', 523, attribute_name )
        ENDIF
 #endif
@@ -4022,5 +3803,5 @@
        IF ( global )  THEN
           nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
-          CALL handle_error( 'get_attribute_real global', 524 )
+          CALL handle_error( 'get_attribute_real global', 524, attribute_name )
 !
 !-- Read attributes referring to a single variable. Therefore, first inquire
@@ -4028,7 +3809,7 @@
        ELSE
           nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
-          CALL handle_error( 'get_attribute_real', 524 )
+          CALL handle_error( 'get_attribute_real', 524, attribute_name )
           nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
-          CALL handle_error( 'get_attribute_real', 524 )
+          CALL handle_error( 'get_attribute_real', 524, attribute_name )
        ENDIF
 #endif
@@ -4063,5 +3844,5 @@
        IF ( global )  THEN
           nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
-          CALL handle_error( 'get_attribute_string global', 525 )
+          CALL handle_error( 'get_attribute_string global', 525, attribute_name )
 !
 !--    Read attributes referring to a single variable. Therefore, first inquire
@@ -4069,8 +3850,8 @@
        ELSE
           nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
-          CALL handle_error( 'get_attribute_string', 525 )
+          CALL handle_error( 'get_attribute_string', 525, attribute_name )
 
           nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
-          CALL handle_error( 'get_attribute_string',525 )
+          CALL handle_error( 'get_attribute_string',525, attribute_name )
 
        ENDIF
@@ -4102,9 +3883,9 @@
 !--    First, inquire dimension ID
        nc_stat = NF90_INQ_DIMID( id, TRIM( variable_name ), id_dim )
-       CALL handle_error( 'get_dimension_length', 526 )
+       CALL handle_error( 'get_dimension_length', 526, variable_name )
 !
 !--    Inquire dimension length
        nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, dum, LEN = dim_len )
-       CALL handle_error( 'get_dimension_length', 526 )
+       CALL handle_error( 'get_dimension_length', 526, variable_name )
 
 #endif
@@ -4133,9 +3914,9 @@
 !--    First, inquire variable ID
        nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
-       CALL handle_error( 'get_variable_1d_int', 527 )
+       CALL handle_error( 'get_variable_1d_int', 527, variable_name )
 !
 !--    Inquire dimension length
        nc_stat = NF90_GET_VAR( id, id_var, var )
-       CALL handle_error( 'get_variable_1d_int', 527 )
+       CALL handle_error( 'get_variable_1d_int', 527, variable_name )
 
 #endif
@@ -4164,9 +3945,9 @@
 !--    First, inquire variable ID
        nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
-       CALL handle_error( 'get_variable_1d_real', 527 )
+       CALL handle_error( 'get_variable_1d_real', 528, variable_name )
 !
 !--    Inquire dimension length
        nc_stat = NF90_GET_VAR( id, id_var, var )
-       CALL handle_error( 'get_variable_1d_real', 527 )
+       CALL handle_error( 'get_variable_1d_real', 528, variable_name )
 
 #endif
@@ -4208,5 +3989,5 @@
                                start = (/ 1,      t     /),                    &
                                count = (/ n_file, 1     /) )
-       CALL handle_error( 'get_variable_pr', 527 )
+       CALL handle_error( 'get_variable_pr', 529, variable_name )
 
 #endif
@@ -4220,5 +4001,5 @@
 !> i.e. each core reads its own domain in slices along x.
 !------------------------------------------------------------------------------!
-    SUBROUTINE get_variable_2d_real( id, variable_name, i, var )
+    SUBROUTINE get_variable_2d_real( id, variable_name, var, is, ie, js, je )
 
        USE indices
@@ -4229,9 +4010,16 @@
        CHARACTER(LEN=*)              ::  variable_name   !< variable name
 
-       INTEGER(iwp), INTENT(IN)      ::  i               !< index along x direction
+       INTEGER(iwp)                  ::  i               !< running index along x direction
+       INTEGER(iwp)                  ::  ie              !< start index for subdomain input along x direction
+       INTEGER(iwp)                  ::  is              !< end index for subdomain input along x direction
        INTEGER(iwp), INTENT(IN)      ::  id              !< file id
        INTEGER(iwp)                  ::  id_var          !< variable id
-
-       REAL(wp), DIMENSION(nys:nyn), INTENT(INOUT) ::  var  !< variable to be read
+       INTEGER(iwp)                  ::  j               !< running index along y direction
+       INTEGER(iwp)                  ::  je              !< start index for subdomain input along y direction
+       INTEGER(iwp)                  ::  js              !< end index for subdomain input along y direction
+       
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  tmp   !< temporary variable to read data from file according
+                                                         !< to its reverse memory access
+       REAL(wp), DIMENSION(:,:), INTENT(INOUT) ::  var   !< variable to be read
 #if defined( __netcdf )
 !
@@ -4239,10 +4027,29 @@
        nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
 !
+!--    Check for collective read-operation and set respective NetCDF flags if
+!--    required. 
+       IF ( collective_read )  THEN
+          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
+       ENDIF
+!
+!--    Allocate temporary variable according to memory access on file. 
+       ALLOCATE( tmp(is:ie,js:je) )
+!
 !--    Get variable
-       nc_stat = NF90_GET_VAR( id, id_var, var(nys:nyn),                       &
-                               start = (/ i+1, nys+1 /),                       &
-                               count = (/ 1, nyn - nys + 1 /) )
-
-       CALL handle_error( 'get_variable_2d_real', 528 )
+       nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
+                               start = (/ is+1,      js+1 /),                  &
+                               count = (/ ie-is + 1, je-js+1 /) )   
+                               
+       CALL handle_error( 'get_variable_2d_real', 530, variable_name )
+!
+!--    Resort data. Please note, dimension subscripts of var all start at 1. 
+       DO  i = is, ie 
+          DO  j = js, je 
+             var(j-js+1,i-is+1) = tmp(i,j)
+          ENDDO
+       ENDDO
+       
+       DEALLOCATE( tmp )
+
 #endif
     END SUBROUTINE get_variable_2d_real
@@ -4254,5 +4061,5 @@
 !> i.e. each core reads its own domain in slices along x.
 !------------------------------------------------------------------------------!
-    SUBROUTINE get_variable_2d_int32( id, variable_name, i, var )
+    SUBROUTINE get_variable_2d_int32( id, variable_name, var, is, ie, js, je )
 
        USE indices
@@ -4263,8 +4070,16 @@
        CHARACTER(LEN=*)              ::  variable_name   !< variable name
 
-       INTEGER(iwp), INTENT(IN)      ::  i               !< index along x direction
+       INTEGER(iwp)                  ::  i               !< running index along x direction
+       INTEGER(iwp)                  ::  ie              !< start index for subdomain input along x direction
+       INTEGER(iwp)                  ::  is              !< end index for subdomain input along x direction
        INTEGER(iwp), INTENT(IN)      ::  id              !< file id
        INTEGER(iwp)                  ::  id_var          !< variable id
-       INTEGER(iwp), DIMENSION(nys:nyn), INTENT(INOUT) ::  var  !< variable to be read
+       INTEGER(iwp)                  ::  j               !< running index along y direction
+       INTEGER(iwp)                  ::  je              !< start index for subdomain input along y direction
+       INTEGER(iwp)                  ::  js              !< end index for subdomain input along y direction
+       
+       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  tmp  !< temporary variable to read data from file according
+                                                            !< to its reverse memory access
+       INTEGER(iwp), DIMENSION(:,:), INTENT(INOUT) ::  var  !< variable to be read
 #if defined( __netcdf )
 !
@@ -4272,10 +4087,29 @@
        nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
 !
+!--    Check for collective read-operation and set respective NetCDF flags if
+!--    required. 
+       IF ( collective_read )  THEN
+          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
+       ENDIF
+!
+!--    Allocate temporary variable according to memory access on file. 
+       ALLOCATE( tmp(is:ie,js:je) )
+!
 !--    Get variable
-       nc_stat = NF90_GET_VAR( id, id_var, var(nys:nyn),                       &
-                               start = (/ i+1, nys+1 /),                       &
-                               count = (/ 1, nyn - nys + 1 /) )
-
-       CALL handle_error( 'get_variable_2d_int32', 529 )
+       nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
+                               start = (/ is+1,      js+1 /),                  &
+                               count = (/ ie-is + 1, je-js+1 /) )    
+                               
+       CALL handle_error( 'get_variable_2d_int32', 531, variable_name )                             
+!
+!--    Resort data. Please note, dimension subscripts of var all start at 1. 
+       DO  i = is, ie 
+          DO  j = js, je 
+             var(j-js+1,i-is+1) = tmp(i,j)
+          ENDDO
+       ENDDO
+       
+       DEALLOCATE( tmp )
+
 #endif
     END SUBROUTINE get_variable_2d_int32
@@ -4287,5 +4121,5 @@
 !> i.e. each core reads its own domain in slices along x.
 !------------------------------------------------------------------------------!
-    SUBROUTINE get_variable_2d_int8( id, variable_name, i, var )
+    SUBROUTINE get_variable_2d_int8( id, variable_name, var, is, ie, js, je )
 
        USE indices
@@ -4296,8 +4130,16 @@
        CHARACTER(LEN=*)              ::  variable_name   !< variable name
 
-       INTEGER(iwp), INTENT(IN)      ::  i               !< index along x direction
+       INTEGER(iwp)                  ::  i               !< running index along x direction
+       INTEGER(iwp)                  ::  ie              !< start index for subdomain input along x direction
+       INTEGER(iwp)                  ::  is              !< end index for subdomain input along x direction
        INTEGER(iwp), INTENT(IN)      ::  id              !< file id
        INTEGER(iwp)                  ::  id_var          !< variable id
-       INTEGER(KIND=1), DIMENSION(nys:nyn), INTENT(INOUT) ::  var  !< variable to be read
+       INTEGER(iwp)                  ::  j               !< running index along y direction
+       INTEGER(iwp)                  ::  je              !< start index for subdomain input along y direction
+       INTEGER(iwp)                  ::  js              !< end index for subdomain input along y direction
+       
+       INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE   ::  tmp  !< temporary variable to read data from file according
+                                                               !< to its reverse memory access
+       INTEGER(KIND=1), DIMENSION(:,:), INTENT(INOUT) ::  var  !< variable to be read
 #if defined( __netcdf )
 !
@@ -4305,12 +4147,32 @@
        nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
 !
+!--    Check for collective read-operation and set respective NetCDF flags if
+!--    required. 
+       IF ( collective_read )  THEN
+          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
+       ENDIF
+!
+!--    Allocate temporary variable according to memory access on file. 
+       ALLOCATE( tmp(is:ie,js:je) )
+!
 !--    Get variable
-       nc_stat = NF90_GET_VAR( id, id_var, var(nys:nyn),                       &
-                               start = (/ i+1, nys+1 /),                       &
-                               count = (/ 1, nyn - nys + 1 /) )
-
-       CALL handle_error( 'get_variable_2d_int8', 530 )
+       nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
+                               start = (/ is+1,      js+1 /),                  &
+                               count = (/ ie-is + 1, je-js+1 /) )   
+                               
+       CALL handle_error( 'get_variable_2d_int8', 532, variable_name )
+!
+!--    Resort data. Please note, dimension subscripts of var all start at 1. 
+       DO  i = is, ie 
+          DO  j = js, je 
+             var(j-js+1,i-is+1) = tmp(i,j)
+          ENDDO
+       ENDDO
+       
+       DEALLOCATE( tmp )
+
 #endif
     END SUBROUTINE get_variable_2d_int8
+
 
 !------------------------------------------------------------------------------!
@@ -4319,5 +4181,6 @@
 !> Reads a 3D 8-bit INTEGER variable from file.
 !------------------------------------------------------------------------------!
-    SUBROUTINE get_variable_3d_int8( id, variable_name, i, j, var )
+    SUBROUTINE get_variable_3d_int8( id, variable_name, var, is, ie, js, je,   &
+                                     ks, ke )
 
        USE indices
@@ -4328,30 +4191,53 @@
        CHARACTER(LEN=*)              ::  variable_name   !< variable name
 
-       INTEGER(iwp), INTENT(IN)      ::  i               !< index along x direction
+       INTEGER(iwp)                  ::  i               !< index along x direction
+       INTEGER(iwp)                  ::  ie              !< start index for subdomain input along x direction
+       INTEGER(iwp)                  ::  is              !< end index for subdomain input along x direction
        INTEGER(iwp), INTENT(IN)      ::  id              !< file id
        INTEGER(iwp)                  ::  id_var          !< variable id
-       INTEGER(iwp), INTENT(IN)      ::  j               !< index along y direction
-       INTEGER(iwp)                  ::  n_file          !< number of data-points along 3rd dimension
-
-       INTEGER(iwp), DIMENSION(1:3)  ::  id_dim
-
-       INTEGER( KIND = 1 ), DIMENSION(nzb:nzt+1), INTENT(INOUT) ::  var  !< variable to be read
+       INTEGER(iwp)                  ::  j               !< index along y direction
+       INTEGER(iwp)                  ::  je              !< start index for subdomain input along y direction
+       INTEGER(iwp)                  ::  js              !< end index for subdomain input along y direction
+       INTEGER(iwp)                  ::  k               !< index along any 3rd dimension
+       INTEGER(iwp)                  ::  ke              !< start index of 3rd dimension
+       INTEGER(iwp)                  ::  ks              !< end index of 3rd dimension
+       
+       INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp  !< temporary variable to read data from file according 
+                                                                 !< to its reverse memory access
+
+       INTEGER(KIND=1), DIMENSION(:,:,:), INTENT(INOUT) ::  var  !< variable to be read
 #if defined( __netcdf )
 
 !
 !--    Inquire variable id
-       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
-!
-!--    Get length of first dimension, required for the count parameter.
-!--    Therefore, first inquired dimension ids
-       nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
-       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(3), LEN = n_file )
+       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )   
+!
+!--    Check for collective read-operation and set respective NetCDF flags if
+!--    required. 
+       IF ( collective_read )  THEN
+          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
+       ENDIF
+!
+!--    Allocate temporary variable according to memory access on file. 
+       ALLOCATE( tmp(is:ie,js:je,ks:ke) )
 !
 !--    Get variable
-       nc_stat = NF90_GET_VAR( id, id_var, var,                  &
-                               start = (/ i+1, j+1, 1 /),                      &
-                               count = (/ 1, 1, n_file /) )
-
-       CALL handle_error( 'get_variable_3d_int8', 531 )
+       nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
+                               start = (/ is+1,    js+1,    ks+1 /),           &
+                               count = (/ ie-is+1, je-js+1, ke-ks+1 /) )                              
+
+       CALL handle_error( 'get_variable_3d_int8', 533, variable_name )                               
+!
+!--    Resort data. Please note, dimension subscripts of var all start at 1. 
+       DO  i = is, ie 
+          DO  j = js, je
+             DO  k = ks, ke
+                var(k-ks+1,j-js+1,i-is+1) = tmp(i,j,k)
+             ENDDO
+          ENDDO
+       ENDDO
+       
+       DEALLOCATE( tmp )
+
 #endif
     END SUBROUTINE get_variable_3d_int8
@@ -4363,5 +4249,6 @@
 !> Reads a 3D float variable from file.
 !------------------------------------------------------------------------------!
-    SUBROUTINE get_variable_3d_real( id, variable_name, i, j, var )
+    SUBROUTINE get_variable_3d_real( id, variable_name, var, is, ie, js, je,   &
+                                     ks, ke )
 
        USE indices
@@ -4372,30 +4259,52 @@
        CHARACTER(LEN=*)              ::  variable_name   !< variable name
 
-       INTEGER(iwp), INTENT(IN)      ::  i               !< index along x direction
+       INTEGER(iwp)                  ::  i               !< index along x direction
+       INTEGER(iwp)                  ::  ie              !< start index for subdomain input along x direction
+       INTEGER(iwp)                  ::  is              !< end index for subdomain input along x direction
        INTEGER(iwp), INTENT(IN)      ::  id              !< file id
        INTEGER(iwp)                  ::  id_var          !< variable id
-       INTEGER(iwp), INTENT(IN)      ::  j               !< index along y direction
-       INTEGER(iwp)                  ::  n3              !< number of data-points along 3rd dimension
-
-       INTEGER(iwp), DIMENSION(3)    ::  id_dim
-
-       REAL(wp), DIMENSION(:), INTENT(INOUT) ::  var     !< variable to be read
+       INTEGER(iwp)                  ::  j               !< index along y direction
+       INTEGER(iwp)                  ::  je              !< start index for subdomain input along y direction
+       INTEGER(iwp)                  ::  js              !< end index for subdomain input along y direction
+       INTEGER(iwp)                  ::  k               !< index along any 3rd dimension
+       INTEGER(iwp)                  ::  ke              !< start index of 3rd dimension
+       INTEGER(iwp)                  ::  ks              !< end index of 3rd dimension
+       
+       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp !< temporary variable to read data from file according 
+                                                         !< to its reverse memory access
+
+       REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) ::  var  !< variable to be read
 #if defined( __netcdf )
 
 !
 !--    Inquire variable id
-       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
-!
-!--    Get length of first dimension, required for the count parameter.
-!--    Therefore, first inquired dimension ids
-       nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
-       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(3), LEN = n3 )
+       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )  
+!
+!--    Check for collective read-operation and set respective NetCDF flags if
+!--    required. 
+       IF ( collective_read )  THEN
+          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
+       ENDIF
+!
+!--    Allocate temporary variable according to memory access on file. 
+       ALLOCATE( tmp(is:ie,js:je,ks:ke) )
 !
 !--    Get variable
-       nc_stat = NF90_GET_VAR( id, id_var, var,                                &
-                               start = (/ i+1, j+1, 1 /),                      &
-                               count = (/ 1, 1, n3 /) )
-
-       CALL handle_error( 'get_variable_3d_real', 532 )
+       nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
+                               start = (/ is+1,    js+1,    ks+1 /),           &
+                               count = (/ ie-is+1, je-js+1, ke-ks+1 /) )   
+                               
+       CALL handle_error( 'get_variable_3d_real', 534, variable_name )
+!
+!--    Resort data. Please note, dimension subscripts of var all start at 1. 
+       DO  i = is, ie 
+          DO  j = js, je
+             DO  k = ks, ke
+                var(k-ks+1,j-js+1,i-is+1) = tmp(i,j,k)
+             ENDDO
+          ENDDO
+       ENDDO
+       
+       DEALLOCATE( tmp )
 
 #endif
@@ -4407,5 +4316,85 @@
 !> Reads a 3D float array from file.
 !------------------------------------------------------------------------------!
-    SUBROUTINE get_variable_3d_real_v( id, variable_name, is, ie, js, je, var )
+!     SUBROUTINE get_variable_3d_real_v( id, variable_name, is, ie, js, je, var )
+! 
+!        USE indices
+!        USE pegrid
+! 
+!        IMPLICIT NONE
+! 
+!        CHARACTER(LEN=*)              ::  variable_name   !< variable name
+! 
+!        INTEGER(iwp), INTENT(IN)      ::  is,ie           !< index range along x direction
+!        INTEGER(iwp), INTENT(IN)      ::  id              !< file id
+!        INTEGER(iwp)                  ::  id_var          !< variable id
+!        INTEGER(iwp), INTENT(IN)      ::  js,je           !< index range along y direction
+!        INTEGER(iwp)                  ::  n3              !< number of data-points along 3rd dimension
+! 
+!        INTEGER(iwp)                  ::  i,j,k
+!        INTEGER(iwp), DIMENSION(3)    ::  id_dim
+! 
+!        REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) ::  var         !< variable to be read
+! #if defined( __netcdf )
+! !
+! !--    Inside the ...static NetCDF files, the array is stored as float.
+! !--    Therefore single precision is sufficiant for the temporary array
+! 
+!        REAL(sp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_var     !< temporary array to read NetCDF data in i,j,k direction
+! 
+! !kk    Please check, if it is time consuming to do the inquire every time
+! !
+! !--    Inquire variable id
+!        nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
+! !
+! !--    Get length of third dimension, required for the count parameter.
+! !--    Therefore, first inquired dimension ids
+!        nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
+!        nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(3), LEN = n3 )
+! 
+! !
+! !--    Check for collective read-operation and set respective NetCDF flags if
+! !--    required. 
+!        IF ( collective_read )  THEN
+!           nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
+!        ENDIF
+! 
+! !
+! !--    Allocate temporary array ro read NetCDF data in i,j,k direction
+! 
+!        ALLOCATE(tmp_var(is:ie,js:je,n3))
+! !
+! !--    Get variable
+! !--    Read complete local 3-D array in oone call
+! 
+!        nc_stat = NF90_GET_VAR( id, id_var, tmp_var,                            &
+!                                start = (/ is+1, js+1, 1 /),                    &
+!                                count = (/ ie-is+1, je-js+1, n3 /) )
+! 
+!        CALL handle_error( 'get_variable_3d_real', 532 )
+! 
+! !
+! !--    Resort data in k,j,i direction
+! 
+!        DO i=is,ie
+!           DO j=js,je
+!              DO K=1,n3
+!                 var (k,j-js+1,i-is+1) = tmp_var(i,j,k)
+!              END DO
+!           END DO
+!        END DO
+! 
+!        DEALLOCATE(tmp_var)
+! 
+! #endif
+!     END SUBROUTINE get_variable_3d_real_v
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Reads a 4D float variable from file. 
+!------------------------------------------------------------------------------!
+    SUBROUTINE get_variable_4d_real( id, variable_name, var, is, ie, js, je,   &
+                                     k1s, k1e, k2s, k2e )
 
        USE indices
@@ -4416,64 +4405,60 @@
        CHARACTER(LEN=*)              ::  variable_name   !< variable name
 
-       INTEGER(iwp), INTENT(IN)      ::  is,ie           !< index range along x direction
+       INTEGER(iwp)                  ::  i               !< index along x direction
+       INTEGER(iwp)                  ::  ie              !< start index for subdomain input along x direction
+       INTEGER(iwp)                  ::  is              !< end index for subdomain input along x direction
        INTEGER(iwp), INTENT(IN)      ::  id              !< file id
        INTEGER(iwp)                  ::  id_var          !< variable id
-       INTEGER(iwp), INTENT(IN)      ::  js,je           !< index range along y direction
-       INTEGER(iwp)                  ::  n3              !< number of data-points along 3rd dimension
-
-       INTEGER(iwp)                  ::  i,j,k
-       INTEGER(iwp), DIMENSION(3)    ::  id_dim
-
-       REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) ::  var         !< variable to be read
+       INTEGER(iwp)                  ::  j               !< index along y direction
+       INTEGER(iwp)                  ::  je              !< start index for subdomain input along y direction
+       INTEGER(iwp)                  ::  js              !< end index for subdomain input along y direction
+       INTEGER(iwp)                  ::  k1              !< index along 3rd direction
+       INTEGER(iwp)                  ::  k1e             !< start index for 3rd dimension
+       INTEGER(iwp)                  ::  k1s             !< end index for 3rd dimension
+       INTEGER(iwp)                  ::  k2              !< index along 4th direction
+       INTEGER(iwp)                  ::  k2e             !< start index for 4th dimension
+       INTEGER(iwp)                  ::  k2s             !< end index for 4th dimension
+
+       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE   ::  tmp  !< temporary variable to read data from file according
+                                                            !< to its reverse memory access
+       REAL(wp), DIMENSION(:,:,:,:), INTENT(INOUT) ::  var  !< variable to be read
 #if defined( __netcdf )
-!
-!--    Inside the ...static NetCDF files, the array is stored as float.
-!--    Therefore single precision is sufficiant for the temporary array
-
-       REAL(sp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_var     !< temporary array to read NetCDF data in i,j,k direction
-
-!kk    Please check, if it is time consuming to do the inquire every time
+
 !
 !--    Inquire variable id
        nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
 !
-!--    Get length of first dimension, required for the count parameter.
-!--    Therefore, first inquired dimension ids
-       nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
-       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(3), LEN = n3 )
-
-       IF(collective_read)  THEN
+!--    Check for collective read-operation and set respective NetCDF flags if
+!--    required. 
+       IF ( collective_read )  THEN
           nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
        ENDIF
-
-!
-!--    Allocate temporary array ro read NetCDF data in i,j,k direction
-
-       ALLOCATE(tmp_var(is:ie,js:je,n3))
+!
+!--    Allocate temporary variable according to memory access on file. 
+       ALLOCATE( tmp(is:ie,js:je,k1s:k1e,k2s:k2e) )
 !
 !--    Get variable
-!--    Read complete local 3-D array in oone call
-
-       nc_stat = NF90_GET_VAR( id, id_var, tmp_var,                              &
-                               start = (/ is+1, js+1, 1 /),                      &
-                               count = (/ ie-is+1, je-js+1, n3 /) )
-
-       CALL handle_error( 'get_variable_3d_real', 532 )
-
-!
-!--    Resort data in k,j,i direction
-
-       DO i=is,ie
-          DO j=js,je
-             DO K=1,n3
-                var (k,j-js+1,i-is+1) = tmp_var(i,j,k)
-             END DO
-          END DO
-       END DO
-
-       DEALLOCATE(tmp_var)
-
+       nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
+                               start = (/ is+1,    js+1,   k1s+1, k2s+1 /),    &
+                               count = (/ ie-is+1, je-js+1,                    &
+                                          k1e-k1s+1, k2e-k2s+1 /) )
+
+       CALL handle_error( 'get_variable_4d_real', 535, variable_name )
+!
+!--    Resort data. Please note, dimension subscripts of var all start at 1. 
+       DO  i = is, ie 
+          DO  j = js, je
+             DO  k1 = k1s, k1e
+                DO  k2 = k2s, k2e
+                   var(k2-k2s+1,k1-k1s+1,j-js+1,i-is+1) = tmp(i,j,k1,k2)
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDDO
+       
+       DEALLOCATE( tmp )
 #endif
-    END SUBROUTINE get_variable_3d_real_v
+    END SUBROUTINE get_variable_4d_real
+
 
 
@@ -4481,9 +4466,14 @@
 ! Description:
 ! ------------
-!> Reads a 4D float variable from file. Note, in constrast to 3D versions,
-!> dimensions are already inquired and passed so that they are known here.
-!------------------------------------------------------------------------------!
-    SUBROUTINE get_variable_4d_real( id, variable_name, i, j, var, n3, n4 )
-
+!> Reads a 3D float variables from dynamic driver, such as time-dependent xy-, 
+!> xz- or yz-boundary data as well as 3D initialization data. Please note, 
+!> the passed arguments are start indices and number of elements in each 
+!> dimension, which is in contrast to the other 3d versions where start- and 
+!> end indices are passed. The different handling of 3D dynamic variables is 
+!> due to its asymmetry for the u- and v component.
+!------------------------------------------------------------------------------!
+    SUBROUTINE get_variable_3d_real_dynamic( id, variable_name, var,           &
+                            i1s, i2s, i3s, count_1, count_2, count_3, dynamic)
+                                
        USE indices
        USE pegrid
@@ -4493,69 +4483,69 @@
        CHARACTER(LEN=*)              ::  variable_name   !< variable name
 
-       INTEGER(iwp), INTENT(IN)      ::  i               !< index along x direction
+       LOGICAL                       ::  dynamic         !< additional flag just used to select correct overloaded routine from interface block 
+       
+       INTEGER(iwp)                  ::  count_1         !< number of elements to be read along 1st dimension (with respect to file)
+       INTEGER(iwp)                  ::  count_2         !< number of elements to be read along 2nd dimension (with respect to file)
+       INTEGER(iwp)                  ::  count_3         !< number of elements to be read along 3rd dimension (with respect to file)
+       INTEGER(iwp)                  ::  i1              !< running index along 1st dimension on file
+       INTEGER(iwp)                  ::  i1s             !< start index for subdomain input along 1st dimension (with respect to file)
+       INTEGER(iwp)                  ::  i2              !< running index along 2nd dimension on file       
+       INTEGER(iwp)                  ::  i2s             !< start index for subdomain input along 2nd dimension (with respect to file)
+       INTEGER(iwp)                  ::  i3              !< running index along 3rd dimension on file 
+       INTEGER(iwp)                  ::  i3s             !< start index of 3rd dimension, in dynamic file this is either time (2D boundary) or z (3D)
        INTEGER(iwp), INTENT(IN)      ::  id              !< file id
        INTEGER(iwp)                  ::  id_var          !< variable id
-       INTEGER(iwp), INTENT(IN)      ::  j               !< index along y direction
-       INTEGER(iwp), INTENT(IN)      ::  n3              !< number of data-points along 3rd dimension
-       INTEGER(iwp), INTENT(IN)      ::  n4              !< number of data-points along 4th dimension
-
-       INTEGER(iwp), DIMENSION(3)    ::  id_dim
-
-       REAL(wp), DIMENSION(:,:), INTENT(INOUT) ::  var     !< variable to be read
+       INTEGER(iwp)                  ::  lb1             !< lower bound of 1st dimension (with respect to file)
+       INTEGER(iwp)                  ::  lb2             !< lower bound of 2nd dimension (with respect to file)
+       INTEGER(iwp)                  ::  lb3             !< lower bound of 3rd dimension (with respect to file)
+       INTEGER(iwp)                  ::  ub1             !< upper bound of 1st dimension (with respect to file)
+       INTEGER(iwp)                  ::  ub2             !< upper bound of 2nd dimension (with respect to file)
+       INTEGER(iwp)                  ::  ub3             !< upper bound of 3rd dimension (with respect to file)
+
+       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp !< temporary variable to read data from file according
+                                                         !< to its reverse memory access
+       
+       REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) ::  var !< input variable
+       
 #if defined( __netcdf )
-
 !
 !--    Inquire variable id
        nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
 !
+!--    Check for collective read-operation and set respective NetCDF flags if
+!--    required. 
+       IF ( collective_read )  THEN
+          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
+       ENDIF    
+!
+!--    Allocate temporary variable according to memory access on file. 
+!--    Therefore, determine dimension bounds of input array. 
+       lb1 = LBOUND(var,3)
+       ub1 = UBOUND(var,3)
+       lb2 = LBOUND(var,2)
+       ub2 = UBOUND(var,2)
+       lb3 = LBOUND(var,1)
+       ub3 = UBOUND(var,1)
+       ALLOCATE( tmp(lb1:ub1,lb2:ub2,lb3:ub3) )
+!
 !--    Get variable
-       nc_stat = NF90_GET_VAR( id, id_var, var,                                &
-                               start = (/ i+1, j+1, 1, 1 /),                   &
-                               count = (/ 1, 1, n3, n4 /) )
-
-       CALL handle_error( 'get_variable_4d_real', 533 )
+       nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
+                               start = (/ i1s,     i2s,     i3s /),            &
+                               count = (/ count_1, count_2, count_3 /) )
+
+       CALL handle_error( 'get_variable_3d_real_dynamic', 536, variable_name )
+!
+!--    Resort data. Please note, dimension subscripts of var all start at 1. 
+       DO  i3 = lb3, ub3
+          DO i2 = lb2, ub2
+             DO  i1 = lb1, ub1
+                var(i3,i2,i1) = tmp(i1,i2,i3)
+             ENDDO
+          ENDDO
+       ENDDO
+       
+       DEALLOCATE( tmp )       
 #endif
-    END SUBROUTINE get_variable_4d_real
-
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Reads a 3D float variable at left, right, north, south and top boundaries.
-!------------------------------------------------------------------------------!
-    SUBROUTINE get_variable_bc( id, variable_name, t_start,                    &
-                                i2_s, count_2, i3_s, count_3,  var )
-
-       USE indices
-       USE pegrid
-
-       IMPLICIT NONE
-
-       CHARACTER(LEN=*)              ::  variable_name   !< variable name
-
-       INTEGER(iwp)                  ::  count_2         !< number of elements in second dimension
-       INTEGER(iwp)                  ::  count_3         !< number of elements in third dimension (usually 1)
-       INTEGER(iwp)                  ::  i2_s            !< start index of second dimension
-       INTEGER(iwp)                  ::  i3_s            !< start index of third dimension
-       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
-       INTEGER(iwp)                  ::  id_var          !< variable id
-       INTEGER(iwp)                  ::  t_start         !< start index at time dimension with respect to netcdf convention
-
-       REAL(wp), DIMENSION(:), INTENT(INOUT) ::  var     !< input variable
-#if defined( __netcdf )
-
-!
-!--    Inquire variable id
-       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
-!
-!--    Get variable
-       nc_stat = NF90_GET_VAR( id, id_var, var,                              &
-                               start = (/ i3_s, i2_s, t_start /),            &
-                               count = (/ count_3, count_2, 1 /) )
-
-       CALL handle_error( 'get_variable_bc', 532 )
-#endif
-    END SUBROUTINE get_variable_bc
+    END SUBROUTINE get_variable_3d_real_dynamic
 
 
@@ -4578,5 +4568,5 @@
 
        nc_stat = NF90_INQUIRE( id, NVARIABLES = num_vars )
-       CALL handle_error( 'inquire_num_variables', 534 )
+       CALL handle_error( 'inquire_num_variables', 537 )
 
 #endif
@@ -4605,9 +4595,9 @@
        ALLOCATE( varids(1:SIZE(var_names)) )
        nc_stat = NF90_INQ_VARIDS( id, NVARS = num_vars, VARIDS = varids )
-       CALL handle_error( 'inquire_variable_names', 535 )
+       CALL handle_error( 'inquire_variable_names', 538 )
 
        DO  i = 1, SIZE(var_names)
           nc_stat = NF90_INQUIRE_VARIABLE( id, varids(i), NAME = var_names(i) )
-          CALL handle_error( 'inquire_variable_names', 535 )
+          CALL handle_error( 'inquire_variable_names', 538 )
        ENDDO
 
@@ -4621,5 +4611,5 @@
 !> Prints out a text message corresponding to the current status.
 !------------------------------------------------------------------------------!
-    SUBROUTINE handle_error( routine_name, errno )
+    SUBROUTINE handle_error( routine_name, errno, name )
 
        USE control_parameters,                                                 &
@@ -4630,12 +4620,20 @@
        CHARACTER(LEN=6) ::  message_identifier
        CHARACTER(LEN=*) ::  routine_name
+       CHARACTER(LEN=*), OPTIONAL ::  name
 
        INTEGER(iwp) ::  errno
 #if defined( __netcdf )
-
+       
        IF ( nc_stat /= NF90_NOERR )  THEN
 
           WRITE( message_identifier, '(''NC'',I4.4)' )  errno
-          message_string = TRIM( NF90_STRERROR( nc_stat ) )
+          
+          IF ( PRESENT( name ) )  THEN
+             message_string = "Problem reading attribute/variable - " //       &
+                              TRIM(name) // ": " //                            &
+                              TRIM( NF90_STRERROR( nc_stat ) )
+          ELSE
+             message_string = TRIM( NF90_STRERROR( nc_stat ) )
+          ENDIF
 
           CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
Index: /palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
===================================================================
--- /palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 3050)
+++ /palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90	(revision 3051)
@@ -25,8 +25,8 @@
 ! -----------------
 ! $Id$
-! Error messages revised
+! Bugfix in calculation of initial Reynolds-stress tensor.
 ! 
 ! 3045 2018-05-28 07:55:41Z Giersch
-! Error message revised
+! Error messages revised
 ! 
 ! 3044 2018-05-25 10:59:41Z gronemeier
@@ -718,8 +718,6 @@
        IF ( a22(k) > 0.0_wp )  THEN
           a22(k) = SQRT( a22(k) )
-          a32(k) = ( r32(k) - a21(k) * a31(k) ) / a22(k)
        ELSE
           a22(k) = 0.0_wp
-          a32(k) = 0.0_wp
        ENDIF
 
