Index: palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90
===================================================================
--- palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90	(revision 2342)
+++ palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90	(revision 2696)
@@ -1,5 +1,5 @@
-!> @file lsf_nudging_mod.f90
-!------------------------------------------------------------------------------!
-! This file is part of PALM.
+!> @file large_scale_forcing_nudging_mod.f90
+!------------------------------------------------------------------------------!
+! This file is part of the PALM model system.
 !
 ! PALM is free software: you can redistribute it and/or modify it under the
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Forcing with larger-scale models implemented (MS)
+! 
+! 2342 2017-08-08 11:00:43Z boeske
 ! fixed check if surface forcing data is available until end of simulation
 ! 
@@ -46,19 +49,25 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  heatflux_input_conversion, pt, pt_init, q, q_init, tend,        &
-               u, u_init, ug, v, v_init, vg, w_subs,                           &
-               waterflux_input_conversion, zu                  
+        ONLY:  dzw, e, heatflux_input_conversion, pt, pt_init, q, q_init, s,   &
+               tend, u, u_init, ug, v, v_init, vg, w, w_subs,                  &
+               waterflux_input_conversion, zu, zw                  
 
     USE control_parameters,                                                    &
-        ONLY:  bc_lr, bc_ns, bc_pt_b, bc_q_b, data_output_pr, dt_3d, end_time, &  
+        ONLY:  bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion,              &
+               data_output_pr, dt_3d, end_time, forcing,                       &
+               force_bound_l, force_bound_n, force_bound_r, force_bound_s,     &  
                humidity, intermediate_timestep_count, ibc_pt_b, ibc_q_b,       &
                large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
-               lsf_exception, message_string, nudging, passive_scalar,         &
+               lsf_exception, message_string, neutral, nudging, passive_scalar,&
                pt_surface, ocean, q_surface, surface_pressure, topography,     &
                use_subsidence_tendencies
 
+    USE grid_variables
+
+    USE pegrid
 
     USE indices,                                                               &
-        ONLY:  ngp_sums_ls, nxl, nxr, nys, nyn, nzb, nz, nzt, wall_flags_0
+        ONLY:  nbgp, ngp_sums_ls, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,     &
+               nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
 
     USE kinds
@@ -70,7 +79,11 @@
         ONLY:  hom, statistic_regions, sums_ls_l, weight_substep
 
+    USE netcdf_data_input_mod,                                                 &
+        ONLY:  force, netcdf_data_input_interpolate
 
     INTEGER(iwp) ::  nlsf = 1000                       !< maximum number of profiles in LSF_DATA (large scale forcing)
     INTEGER(iwp) ::  ntnudge = 1000                    !< maximum number of profiles in NUDGING_DATA (nudging)
+
+    REAL(wp) ::  d_area_t
 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ptnudge     !< vertical profile of pot. temperature interpolated to vertical grid (nudging)
@@ -107,9 +120,11 @@
            lsf_nudging_check_parameters, nudge_init,                           &
            lsf_nudging_check_data_output_pr, lsf_nudging_header,               &
-           calc_tnudge, nudge, nudge_ref
+           calc_tnudge, nudge, nudge_ref, forcing_bc_mass_conservation,        &
+           forcing_bc
 !
 !-- Public variables
     PUBLIC qsws_surf, shf_surf, td_lsa_lpt, td_lsa_q, td_sub_lpt,              &
-           td_sub_q, time_vert
+           td_sub_q, time_vert, force
+
 
     INTERFACE ls_advec
@@ -125,4 +140,523 @@
  CONTAINS
 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> @todo Missing subroutine description.
+!------------------------------------------------------------------------------!
+    SUBROUTINE forcing_bc_mass_conservation
+
+       USE control_parameters,                                                 &
+           ONLY:  volume_flow
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i !<
+       INTEGER(iwp) ::  j !<
+       INTEGER(iwp) ::  k !<
+
+       REAL(wp) ::  w_correct !< 
+       REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !<
+
+       volume_flow   = 0.0_wp
+       volume_flow_l = 0.0_wp
+
+       d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
+
+       IF ( force_bound_l )  THEN
+          i = nxl
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy   &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,j,i), 1 ) )
+             ENDDO
+          ENDDO
+       ENDIF
+       IF ( force_bound_r )  THEN
+          i = nxr+1
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+                volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy   &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,j,i), 1 ) )
+             ENDDO
+          ENDDO
+       ENDIF
+       IF ( force_bound_s )  THEN
+          j = nys
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt
+                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx   &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,j,i), 2 ) )
+             ENDDO
+          ENDDO
+       ENDIF
+       IF ( force_bound_n )  THEN
+          j = nyn+1
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt
+                volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx   &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,j,i), 2 ) )
+             ENDDO
+          ENDDO
+       ENDIF
+!
+!--       Top boundary
+       k = nzt
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
+          ENDDO
+       ENDDO
+
+#if defined( __parallel )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM,      &
+                           comm2d, ierr )
+#else
+       volume_flow = volume_flow_l
+#endif
+
+       w_correct = SUM( volume_flow ) * d_area_t
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzt, nzt + 1
+                w(k,j,i) = w(k,j,i) + w_correct
+             ENDDO
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE forcing_bc_mass_conservation
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> @todo Missing subroutine description.
+!------------------------------------------------------------------------------!
+    SUBROUTINE forcing_bc
+
+       USE control_parameters,                                                 &
+           ONLY:  force_bound_l, force_bound_n, force_bound_r, force_bound_s,  &
+                  humidity, neutral, passive_scalar, simulated_time
+
+       USE netcdf_data_input_mod,                                              &
+           ONLY:  force                              
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i !< running index x-direction
+       INTEGER(iwp) ::  j !< running index y-direction
+       INTEGER(iwp) ::  k !< running index z-direction
+       INTEGER(iwp) ::  t !< running index for time levels
+
+       REAL(wp) ::  ddt_lsf !< inverse value of time resolution of forcing data
+       REAL(wp) ::  t_ref   !< time past since last reference step
+
+!
+!--    If required, interpolate and/or extrapolate data vertically. This is 
+!--    required as Inifor outputs only equidistant vertical data.
+       IF ( ANY( zu(1:nzt+1) /= force%zu_atmos(1:force%nzu) ) )  THEN
+          IF ( .NOT. force%interpolated )  THEN
+
+             DO  t = 0, 1
+                IF ( force_bound_l )  THEN
+                   CALL netcdf_data_input_interpolate( force%u_left(t,:,:),    &
+                                                       zu(nzb+1:nzt+1),        &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%v_left(t,:,:),    &
+                                                       zu(nzb+1:nzt+1),        &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%w_left(t,:,:),    &
+                                                       zw(nzb+1:nzt+1),        &
+                                                       force%zw_atmos )
+                   IF ( .NOT. neutral )                                        &
+                      CALL netcdf_data_input_interpolate( force%pt_left(t,:,:),&
+                                                          zu(nzb+1:nzt+1),     &
+                                                          force%zu_atmos )
+                   IF ( humidity )                                             &
+                      CALL netcdf_data_input_interpolate( force%q_left(t,:,:), &
+                                                          zu(nzb+1:nzt+1),     &
+                                                          force%zu_atmos )
+                ENDIF
+                IF ( force_bound_r )  THEN
+                   CALL netcdf_data_input_interpolate( force%u_right(t,:,:),   &
+                                                       zu(nzb+1:nzt+1),        &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%v_right(t,:,:),   &
+                                                       zu(nzb+1:nzt+1),        &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%w_right(t,:,:),   &
+                                                       zw(nzb+1:nzt+1),        &
+                                                       force%zw_atmos )
+                   IF ( .NOT. neutral )                                        &
+                      CALL netcdf_data_input_interpolate( force%pt_right(t,:,:),&
+                                                          zu(nzb+1:nzt+1),     &
+                                                          force%zu_atmos )
+                   IF ( humidity )                                             &
+                      CALL netcdf_data_input_interpolate( force%q_right(t,:,:),&
+                                                          zu(nzb+1:nzt+1),     &
+                                                          force%zu_atmos )
+                ENDIF
+                IF ( force_bound_n )  THEN
+                   CALL netcdf_data_input_interpolate( force%u_north(t,:,:),   &
+                                                       zu(nzb+1:nzt+1),        &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%v_north(t,:,:),   &
+                                                       zu(nzb+1:nzt+1),        &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%w_north(t,:,:),   &
+                                                       zw(nzb+1:nzt+1),        &
+                                                       force%zw_atmos )
+                   IF ( .NOT. neutral )                                        &
+                      CALL netcdf_data_input_interpolate( force%pt_north(t,:,:),&
+                                                          zu(nzb+1:nzt+1),     &
+                                                          force%zu_atmos )
+                   IF ( humidity )                                             &
+                      CALL netcdf_data_input_interpolate( force%q_north(t,:,:),&
+                                                          zu(nzb+1:nzt+1),     &
+                                                          force%zu_atmos )
+                ENDIF
+                IF ( force_bound_s )  THEN
+                   CALL netcdf_data_input_interpolate( force%u_south(t,:,:),   &
+                                                       zu(nzb+1:nzt+1),        &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%v_south(t,:,:),   &
+                                                       zu(nzb+1:nzt+1),        &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%w_south(t,:,:),   &
+                                                       zw(nzb+1:nzt+1),        &
+                                                       force%zw_atmos )
+                   IF ( .NOT. neutral )                                        &
+                      CALL netcdf_data_input_interpolate( force%pt_south(t,:,:),&
+                                                          zu(nzb+1:nzt+1),     &
+                                                          force%zu_atmos )
+                   IF ( humidity )                                             &
+                      CALL netcdf_data_input_interpolate( force%q_south(t,:,:),&
+                                                          zu(nzb+1:nzt+1),     &
+                                                          force%zu_atmos )
+                ENDIF
+             ENDDO
+!
+!--          Note, no interpolation of top boundary. Just use initial value. 
+!--          No physical meaningful extrapolation possible if only one layer is
+!--          given. 
+
+             force%interpolated = .TRUE. 
+          ENDIF
+       ENDIF
+      
+!
+!--    Calculate time interval of forcing data       
+       ddt_lsf = 1.0_wp / ( force%time(force%tind_p) - force%time(force%tind) )
+!
+!--    Calculate reziproke time past since last reference step. Please note, 
+!--    as simulated time is still not updated, the actual time here is 
+!--    simulated time + dt_3d
+       t_ref = simulated_time + dt_3d - force%time(force%tind)
+
+       IF ( force_bound_l )  THEN
+
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt+1
+                u(k,j,nxlg:nxl) = force%u_left(0,k,j) + ddt_lsf * t_ref *      &
+                         ( force%u_left(1,k,j) - force%u_left(0,k,j) ) *       &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,j,nxlg:nxl), 1 ) )
+             ENDDO
+          ENDDO
+
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+                w(k,j,nxlg:nxl-1) = force%w_left(0,k,j) + ddt_lsf * t_ref *    &
+                         ( force%w_left(1,k,j) - force%w_left(0,k,j) ) *       &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,j,nxlg:nxl-1), 3 ) )
+             ENDDO
+          ENDDO
+
+          DO  j = nysv, nyn
+             DO  k = nzb+1, nzt+1
+                v(k,j,nxlg:nxl-1) = force%v_left(0,k,j) + ddt_lsf * t_ref *    &
+                         ( force%v_left(1,k,j) - force%v_left(0,k,j) ) *       &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,j,nxlg:nxl-1), 2 ) )
+             ENDDO
+          ENDDO
+
+          IF ( .NOT. neutral )  THEN
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt+1
+                   pt(k,j,nxlg:nxl-1) = force%pt_left(0,k,j) + ddt_lsf *       &
+                                                                  t_ref   *    &
+                       ( force%pt_left(1,k,j) - force%pt_left(0,k,j) )
+ 
+                ENDDO
+             ENDDO
+          ENDIF
+
+          IF ( humidity )  THEN
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt+1
+                   q(k,j,nxlg:nxl-1) = force%q_left(0,k,j) + ddt_lsf   *       &
+                                                                  t_ref   *    &
+                       ( force%q_left(1,k,j) - force%q_left(0,k,j) )
+ 
+                ENDDO
+             ENDDO
+          ENDIF
+
+       ENDIF
+
+       IF ( force_bound_r )  THEN
+
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt+1
+                u(k,j,nxr+1:nxrg) = force%u_right(0,k,j) + ddt_lsf * t_ref *   &
+                        ( force%u_right(1,k,j) - force%u_right(0,k,j) ) *      &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,j,nxr+1:nxrg), 1 ) )
+
+             ENDDO
+          ENDDO
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+                w(k,j,nxr+1:nxrg) = force%w_right(0,k,j) + ddt_lsf * t_ref *   &
+                        ( force%w_right(1,k,j) - force%w_right(0,k,j) ) *      &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,j,nxr+1:nxrg), 3 ) )
+             ENDDO
+          ENDDO
+
+          DO  j = nysv, nyn
+             DO  k = nzb+1, nzt+1
+                v(k,j,nxr+1:nxrg) = force%v_right(0,k,j) + ddt_lsf * t_ref *   &
+                        ( force%v_right(1,k,j) - force%v_right(0,k,j) ) *      &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,j,nxr+1:nxrg), 2 ) )
+             ENDDO
+          ENDDO
+
+          IF ( .NOT. neutral )  THEN
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt+1
+                   pt(k,j,nxr+1:nxrg) = force%pt_right(0,k,j) + ddt_lsf *      &
+                                                                  t_ref   *    &
+                     ( force%pt_right(1,k,j) - force%pt_right(0,k,j) )
+ 
+                ENDDO
+             ENDDO
+          ENDIF
+
+          IF ( humidity )  THEN
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt+1
+                   q(k,j,nxr+1:nxrg) = force%q_right(0,k,j) + ddt_lsf   *      &
+                                                                  t_ref   *    &
+                       ( force%q_right(1,k,j) - force%q_right(0,k,j) )
+ 
+                ENDDO
+             ENDDO
+          ENDIF
+
+       ENDIF
+
+       IF ( force_bound_s )  THEN
+
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt+1
+                v(k,nysg:nys,i)   = force%v_south(0,k,i) + ddt_lsf * t_ref *   &
+                        ( force%v_south(1,k,i) - force%v_south(0,k,i) ) *      &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,nysg:nys,i), 2 ) )
+             ENDDO
+          ENDDO
+
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt
+                w(k,nysg:nys-1,i) = force%w_south(0,k,i) + ddt_lsf * t_ref *   &
+                        ( force%w_south(1,k,i) - force%w_south(0,k,i) ) *      &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,nysg:nys-1,i), 3 ) )
+             ENDDO
+          ENDDO
+
+          DO  i = nxlu, nxr
+             DO  k = nzb+1, nzt+1
+                u(k,nysg:nys-1,i) = force%u_south(0,k,i) + ddt_lsf * t_ref *   &
+                        ( force%u_south(1,k,i) - force%u_south(0,k,i) ) *      &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,nysg:nys-1,i), 1 ) )
+             ENDDO
+          ENDDO
+
+          IF ( .NOT. neutral )  THEN
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt+1
+                   pt(k,nysg:nys-1,i) = force%pt_south(0,k,i) + ddt_lsf *      &
+                                                                  t_ref   *    &
+                     ( force%pt_south(1,k,i) - force%pt_south(0,k,i) )
+ 
+                ENDDO
+             ENDDO
+          ENDIF
+
+          IF ( humidity )  THEN
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt+1
+                   q(k,nysg:nys-1,i) = force%q_south(0,k,i) + ddt_lsf   *      &
+                                                                  t_ref   *    &
+                       ( force%q_south(1,k,i) - force%q_south(0,k,i) )
+ 
+                ENDDO
+             ENDDO
+          ENDIF
+
+       ENDIF
+
+       IF ( force_bound_n )  THEN
+
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt+1
+                v(k,nyn+1:nyng,i)   = force%v_north(0,k,i) + ddt_lsf * t_ref * &
+                        ( force%v_north(1,k,i) - force%v_north(0,k,i) ) *      &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 2 ) )
+             ENDDO
+          ENDDO
+          DO  i = nxl, nxr
+             DO  k = nzb+1, nzt
+                w(k,nyn+1:nyng,i) = force%w_north(0,k,i) + ddt_lsf * t_ref *   &
+                        ( force%w_north(1,k,i) - force%w_north(0,k,i) ) *      &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 3 ) )
+             ENDDO
+          ENDDO
+
+          DO  i = nxlu, nxr
+             DO  k = nzb+1, nzt+1
+                u(k,nyn+1:nyng,i) = force%u_north(0,k,i) + ddt_lsf * t_ref *   &
+                        ( force%u_north(1,k,i) - force%u_north(0,k,i) ) *      &
+                           MERGE( 1.0_wp, 0.0_wp,                              &
+                                  BTEST( wall_flags_0(k,nyn+1:nyng,i), 1 ) )
+
+             ENDDO
+          ENDDO
+
+          IF ( .NOT. neutral )  THEN
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt+1
+                   pt(k,nyn+1:nyng,i) = force%pt_north(0,k,i) + ddt_lsf *      &
+                                                                  t_ref   *    &
+                     ( force%pt_north(1,k,i) - force%pt_north(0,k,i) )
+ 
+                ENDDO
+             ENDDO
+          ENDIF
+
+          IF ( humidity )  THEN
+             DO  i = nxl, nxr
+                DO  k = nzb+1, nzt+1
+                   q(k,nyn+1:nyng,i) = force%q_north(0,k,i) + ddt_lsf   *      &
+                                                                  t_ref   *    &
+                       ( force%q_north(1,k,i) - force%q_north(0,k,i) )
+ 
+                ENDDO
+             ENDDO
+          ENDIF
+
+       ENDIF
+!
+!--    Top boundary. 
+!--    Please note, only map Inifor data on model top in case the numeric is
+!--    identical to the Inifor grid. At the top boundary an extrapolation is 
+!--    not possible. 
+       IF ( ANY( zu(1:nzt+1) /= force%zu_atmos(1:force%nzu) ) )  THEN
+          DO  i = nxlu, nxr
+             DO  j = nys, nyn
+                u(nzt+1,j,i) = force%u_top(0,j,i) + ddt_lsf * t_ref *          &
+                           ( force%u_top(1,j,i) - force%u_top(0,j,i) ) *       &
+                              MERGE( 1.0_wp, 0.0_wp,                           &
+                                     BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
+             ENDDO
+          ENDDO
+
+          DO  i = nxl, nxr
+             DO  j = nysv, nyn
+                v(nzt+1,j,i) = force%v_top(0,j,i) + ddt_lsf * t_ref *          &
+                           ( force%v_top(1,j,i) - force%v_top(0,j,i) ) *       &
+                              MERGE( 1.0_wp, 0.0_wp,                           &
+                                     BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
+             ENDDO
+          ENDDO
+
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                w(nzt:nzt+1,j,i) = force%w_top(0,j,i) + ddt_lsf * t_ref *      &
+                           ( force%w_top(1,j,i) - force%w_top(0,j,i) ) *       &
+                              MERGE( 1.0_wp, 0.0_wp,                           &
+                                     BTEST( wall_flags_0(nzt:nzt+1,j,i), 3 ) )
+             ENDDO
+          ENDDO
+
+
+          IF ( .NOT. neutral )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   pt(nzt+1,j,i) = force%pt_top(0,j,i) + ddt_lsf * t_ref *     &
+                           ( force%pt_top(1,j,i) - force%pt_top(0,j,i) )
+                ENDDO
+             ENDDO
+          ENDIF
+
+          IF ( humidity )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   q(nzt+1,j,i) = force%q_top(0,j,i) + ddt_lsf * t_ref *       &
+                           ( force%q_top(1,j,i) - force%q_top(0,j,i) )
+                ENDDO
+             ENDDO
+          ENDIF
+       ENDIF
+!
+!--    Moreover, set Neumann boundary condition for subgrid-scale TKE and 
+!--    passive scalar
+       IF ( .NOT. constant_diffusion )  THEN
+          IF (  force_bound_l )  e(:,:,nxl-1) = e(:,:,nxl)
+          IF (  force_bound_r )  e(:,:,nxr+1) = e(:,:,nxr)
+          IF (  force_bound_s )  e(:,nys-1,:) = e(:,nys,:)
+          IF (  force_bound_n )  e(:,nyn+1,:) = e(:,nyn,:)
+          e(nzt+1,:,:) = e(nzt,:,:)
+       ENDIF
+       IF ( passive_scalar )  THEN
+          IF (  force_bound_l )  s(:,:,nxl-1) = s(:,:,nxl)
+          IF (  force_bound_r )  s(:,:,nxr+1) = s(:,:,nxr)
+          IF (  force_bound_s )  s(:,nys-1,:) = s(:,nys,:)
+          IF (  force_bound_n )  s(:,nyn+1,:) = s(:,nyn,:)
+       ENDIF
+
+
+       CALL exchange_horiz( u, nbgp )
+       CALL exchange_horiz( v, nbgp )
+       CALL exchange_horiz( w, nbgp )
+       IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
+       IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
+
+!
+!--    Set surface pressure. Please note, time-dependent surface
+!--    pressure would require changes in anelastic approximation and 
+!--    treatment of fluxes. 
+!--    For the moment, comment this out! 
+!      surface_pressure = force%surface_pressure(force%tind) +                 &
+!                                                      ddt_lsf * t_ref *       &
+!                                    ( force%surface_pressure(force%tind_p)    &
+!                                    - force%surface_pressure(force%tind) )
+
+    END SUBROUTINE forcing_bc
 
 !------------------------------------------------------------------------------!
@@ -403,4 +937,7 @@
     SUBROUTINE lsf_init
 
+       USE netcdf_data_input_mod,                                              &
+           ONLY:  netcdf_data_input_lsf  
+
        IMPLICIT NONE
 
@@ -411,5 +948,6 @@
        INTEGER(iwp) ::  finput = 90   !<
        INTEGER(iwp) ::  k             !<
-       INTEGER(iwp) ::  nt             !<
+       INTEGER(iwp) ::  nt            !<
+       INTEGER(iwp) ::  t             !< running index for time levels
 
        REAL(wp) ::  fac               !<
@@ -432,105 +970,224 @@
        REAL(wp) ::  r_dummy           !<
 
-       ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),              &
-                 qsws_surf(0:nlsf), shf_surf(0:nlsf),                          &
-                 td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),     &
-                 td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),     &
-                 time_vert(0:nlsf), time_surf(0:nlsf),                         &
-                 ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),         &
-                 wsubs_vert(nzb:nzt+1,0:nlsf) )
-
-       p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
-       shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
-       td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
-       time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
-       wsubs_vert = 0.0_wp
-
-!
-!--    Array for storing large scale forcing and nudging tendencies at each 
-!--    timestep for data output
-       ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
-       sums_ls_l = 0.0_wp
-
-       ngp_sums_ls = (nz+2)*6
-
-       OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
-              FORM='FORMATTED', IOSTAT=ierrn )
-
-       IF ( ierrn /= 0 )  THEN
-          message_string = 'file LSF_DATA does not exist'
-          CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
-       ENDIF
-
-       ierrn = 0
-!
-!--    First three lines of LSF_DATA contain header
-       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
-       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
-       READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
-
-       IF ( ierrn /= 0 )  THEN
-          message_string = 'errors in file LSF_DATA'
-          CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
-       ENDIF
-
-!
-!--    Surface values are read in
-       nt     = 0
-       ierrn = 0
-
-       DO WHILE ( time_surf(nt) < end_time )
-          nt = nt + 1
-          READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),      &
-                                             qsws_surf(nt), pt_surf(nt),       &
-                                             q_surf(nt), p_surf(nt)
+       IF ( forcing )  THEN
+!
+!--       Allocate arrays for reading boundary values. Arrays will incorporate 2 
+!--       time levels in order to interpolate in between. 
+          IF ( force_bound_l )  THEN
+             ALLOCATE( force%u_left(0:1,nzb+1:nzt+1,nys:nyn)  )
+             ALLOCATE( force%v_left(0:1,nzb+1:nzt+1,nysv:nyn) )
+             ALLOCATE( force%w_left(0:1,nzb+1:nzt,nys:nyn)    )
+             IF ( humidity )      ALLOCATE( force%q_left(0:1,nzb+1:nzt+1,nys:nyn)  )
+             IF ( .NOT. neutral ) ALLOCATE( force%pt_left(0:1,nzb+1:nzt+1,nys:nyn) )
+          ENDIF
+          IF ( force_bound_r )  THEN
+             ALLOCATE( force%u_right(0:1,nzb+1:nzt+1,nys:nyn)  )
+             ALLOCATE( force%v_right(0:1,nzb+1:nzt+1,nysv:nyn) )
+             ALLOCATE( force%w_right(0:1,nzb+1:nzt,nys:nyn)    )
+             IF ( humidity )      ALLOCATE( force%q_right(0:1,nzb+1:nzt+1,nys:nyn)  )
+             IF ( .NOT. neutral ) ALLOCATE( force%pt_right(0:1,nzb+1:nzt+1,nys:nyn) )
+          ENDIF
+          IF ( force_bound_n )  THEN
+             ALLOCATE( force%u_north(0:1,nzb+1:nzt+1,nxlu:nxr) )
+             ALLOCATE( force%v_north(0:1,nzb+1:nzt+1,nxl:nxr)  )
+             ALLOCATE( force%w_north(0:1,nzb+1:nzt,nxl:nxr)    )
+             IF ( humidity )      ALLOCATE( force%q_north(0:1,nzb+1:nzt+1,nxl:nxr)  )
+             IF ( .NOT. neutral ) ALLOCATE( force%pt_north(0:1,nzb+1:nzt+1,nxl:nxr) )
+          ENDIF
+          IF ( force_bound_s )  THEN
+             ALLOCATE( force%u_south(0:1,nzb+1:nzt+1,nxlu:nxr) )
+             ALLOCATE( force%v_south(0:1,nzb+1:nzt+1,nxl:nxr)  )
+             ALLOCATE( force%w_south(0:1,nzb+1:nzt,nxl:nxr)    )
+             IF ( humidity )      ALLOCATE( force%q_south(0:1,nzb+1:nzt+1,nxl:nxr)  )
+             IF ( .NOT. neutral ) ALLOCATE( force%pt_south(0:1,nzb+1:nzt+1,nxl:nxr) )
+          ENDIF
+          
+          ALLOCATE( force%u_top(0:1,nys:nyn,nxlu:nxr) )
+          ALLOCATE( force%v_top(0:1,nysv:nyn,nxl:nxr) )
+          ALLOCATE( force%w_top(0:1,nys:nyn,nxl:nxr)  )
+          IF ( humidity )      ALLOCATE( force%q_top(0:1,nys:nyn,nxl:nxr)  )
+          IF ( .NOT. neutral ) ALLOCATE( force%pt_top(0:1,nys:nyn,nxl:nxr) )
+
+!
+!--       Initial call of input. Time array, initial 3D data of u, v, w, 
+!--       potential temperature, as well as mixing ratio, will be read. 
+!--       Moreover, data at lateral and top boundary will be read.  
+          CALL netcdf_data_input_lsf
+!
+!--       Please note, at the moment INIFOR assumes only an equidistant vertical
+!--       grid. In case of vertical grid stretching, input of inital 3D data 
+!--       need to be inter- and/or extrapolated. 
+!--       Therefore, check if zw grid on file is identical to numeric zw grid.   
+          IF ( ANY( zu(1:nzt+1) /= force%zu_atmos(1:force%nzu) ) )  THEN
+!
+!--          Also data at the boundaries need to be inter/extrapolated at both
+!--          time levels
+             DO  t = 0, 1
+                IF ( force_bound_l )  THEN
+                   CALL netcdf_data_input_interpolate( force%u_left(t,:,:),    &
+                                                       zu(1:nzt+1),            &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%v_left(t,:,:),    &
+                                                       zu(1:nzt+1),            &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%w_left(t,:,:),    &
+                                                       zw(1:nzt+1),            &
+                                                       force%zw_atmos )
+                   IF ( .NOT. neutral )                                        &
+                      CALL netcdf_data_input_interpolate( force%pt_left(t,:,:),&
+                                                          zu(1:nzt+1),         &
+                                                          force%zu_atmos )
+                   IF ( humidity )                                             &
+                      CALL netcdf_data_input_interpolate( force%q_left(t,:,:), &
+                                                          zu(1:nzt+1),         &
+                                                          force%zu_atmos )
+                ENDIF
+                IF ( force_bound_r )  THEN
+                   CALL netcdf_data_input_interpolate( force%u_right(t,:,:),   &
+                                                       zu(1:nzt+1),            &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%v_right(t,:,:),   &
+                                                       zu(1:nzt+1),            &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%w_right(t,:,:),   &
+                                                       zw(1:nzt+1),            &
+                                                       force%zw_atmos )
+                   IF ( .NOT. neutral )                                        &
+                      CALL netcdf_data_input_interpolate( force%pt_right(t,:,:),&
+                                                          zu(1:nzt+1),          &
+                                                          force%zu_atmos )
+                   IF ( humidity )                                             &
+                      CALL netcdf_data_input_interpolate( force%q_right(t,:,:),&
+                                                          zu(1:nzt+1),         &
+                                                          force%zu_atmos )
+                ENDIF
+                IF ( force_bound_n )  THEN
+                   CALL netcdf_data_input_interpolate( force%u_north(t,:,:),   &
+                                                       zu(1:nzt+1),            &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%v_north(t,:,:),   &
+                                                       zu(1:nzt+1),            &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%w_north(t,:,:),   &
+                                                       zw(1:nzt+1),            &
+                                                       force%zw_atmos )
+                   IF ( .NOT. neutral )                                        &
+                      CALL netcdf_data_input_interpolate( force%pt_north(t,:,:),&
+                                                          zu(1:nzt+1),          &
+                                                          force%zu_atmos )
+                   IF ( humidity )                                             &
+                      CALL netcdf_data_input_interpolate( force%q_north(t,:,:),&
+                                                          zu(1:nzt+1),         &
+                                                          force%zu_atmos )
+                ENDIF
+                IF ( force_bound_s )  THEN
+                   CALL netcdf_data_input_interpolate( force%u_south(t,:,:),   &
+                                                       zu(1:nzt+1),            &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%v_south(t,:,:),   &
+                                                       zu(1:nzt+1),            &
+                                                       force%zu_atmos )
+                   CALL netcdf_data_input_interpolate( force%w_south(t,:,:),   &
+                                                       zw(1:nzt+1),            &
+                                                       force%zw_atmos )
+                   IF ( .NOT. neutral )                                        &
+                      CALL netcdf_data_input_interpolate( force%pt_south(t,:,:),&
+                                                          zu(1:nzt+1),          &
+                                                          force%zu_atmos )
+                   IF ( humidity )                                             &
+                      CALL netcdf_data_input_interpolate( force%q_south(t,:,:),&
+                                                          zu(1:nzt+1),         &
+                                                          force%zu_atmos )
+                ENDIF
+             ENDDO
+          ENDIF
+
+!
+!--       Exchange ghost points
+          CALL exchange_horiz( u, nbgp )
+          CALL exchange_horiz( v, nbgp )
+          CALL exchange_horiz( w, nbgp )
+          IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
+          IF ( humidity      )  CALL exchange_horiz( q,  nbgp )
+!
+!--       At lateral boundaries, set also initial boundary conditions
+          IF ( force_bound_l )  THEN
+             u(:,:,nxl)   = u(:,:,nxlu)
+             v(:,:,nxl-1) = v(:,:,nxl)
+             w(:,:,nxl-1) = w(:,:,nxl)
+             IF ( .NOT. neutral )  pt(:,:,nxl-1) = pt(:,:,nxl)
+             IF ( humidity      )  q(:,:,nxl-1)  = q(:,:,nxl)
+          ENDIF
+          IF ( force_bound_r )  THEN
+             u(:,:,nxr+1) = u(:,:,nxr)
+             v(:,:,nxr+1) = v(:,:,nxr)
+             w(:,:,nxr+1) = w(:,:,nxr)
+             IF ( .NOT. neutral )  pt(:,:,nxr+1) = pt(:,:,nxr)
+             IF ( humidity      )  q(:,:,nxr+1)  = q(:,:,nxr)
+          ENDIF
+          IF ( force_bound_s )  THEN
+             u(:,nys-1,:) = u(:,nys,:)
+             v(:,nys,:)   = v(:,nysv,:)
+             w(:,nys-1,:) = w(:,nys,:)
+             IF ( .NOT. neutral )  pt(:,nys-1,:) = pt(:,nys,:)
+             IF ( humidity      )  q(:,nys-1,:)  = q(:,nys,:)
+          ENDIF
+          IF ( force_bound_n )  THEN
+             u(:,nyn+1,:) = u(:,nyn,:)
+             v(:,nyn+1,:) = v(:,nyn,:)
+             w(:,nyn+1,:) = w(:,nyn,:)
+             IF ( .NOT. neutral )  pt(:,nyn+1,:) = pt(:,nyn,:)
+             IF ( humidity      )  q(:,nyn+1,:)  = q(:,nyn,:)
+          ENDIF
+
+!
+!--       After 3D data is initialized, ensure mass conservation
+          CALL forcing_bc_mass_conservation
+!
+!--       Initialize surface pressure. Please note, time-dependent surface
+!--       pressure would require changes in anelastic approximation and 
+!--       treatment of fluxes. 
+!--       For the moment, comment this out! 
+!         surface_pressure = force%surface_pressure(0)
+
+       ELSE
+
+          ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf),           &
+                    qsws_surf(0:nlsf), shf_surf(0:nlsf),                       &
+                    td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf),  &
+                    td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf),  &
+                    time_vert(0:nlsf), time_surf(0:nlsf),                      &
+                    ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf),      &
+                    wsubs_vert(nzb:nzt+1,0:nlsf) )
+
+          p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
+          shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
+          td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
+          time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
+          wsubs_vert = 0.0_wp
+
+!
+!--       Array for storing large scale forcing and nudging tendencies at each 
+!--       timestep for data output
+          ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
+          sums_ls_l = 0.0_wp
+
+          ngp_sums_ls = (nz+2)*6
+
+          OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
+                 FORM='FORMATTED', IOSTAT=ierrn )
 
           IF ( ierrn /= 0 )  THEN
-            WRITE ( message_string, * ) 'No time dependent surface variables ',&
-                              'in&LSF_DATA for end of run found'
-
-             CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
-          ENDIF
-       ENDDO
-
-       IF ( time_surf(1) > end_time )  THEN
-          WRITE ( message_string, * ) 'Time dependent surface variables in ',  &
-                                      '&LSF_DATA set in after end of ' ,      &
-                                      'simulation - lsf_surf is set to FALSE'
-          CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
-          lsf_surf = .FALSE.
-       ENDIF
-
-!
-!--    Go to the end of the list with surface variables
-       DO WHILE ( ierrn == 0 )
-          READ ( finput, *, IOSTAT = ierrn ) r_dummy
-       ENDDO
-
-!
-!--    Profiles of ug, vg and w_subs are read in (large scale forcing)
-
-       nt = 0
-       DO WHILE ( time_vert(nt) < end_time )
-          nt = nt + 1
-          hash = "#"
-          ierrn = 1 ! not zero
-!
-!--       Search for the next line consisting of "# time", 
-!--       from there onwards the profiles will be read
-          DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
-             READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
-             IF ( ierrn < 0 )  THEN 
-                WRITE( message_string, * ) 'No time dependent vertical profiles',&
-                                 ' in&LSF_DATA for end of run found'
-                CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
-             ENDIF
-          ENDDO
-
-          IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
-
-          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,  &
-                                           lowwsubs_vert, low_td_lsa_lpt,      &
-                                           low_td_lsa_q, low_td_sub_lpt,       &
-                                           low_td_sub_q
+             message_string = 'file LSF_DATA does not exist'
+             CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          ierrn = 0
+!
+!--       First three lines of LSF_DATA contain header
+          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
+          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
+          READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
+
           IF ( ierrn /= 0 )  THEN
              message_string = 'errors in file LSF_DATA'
@@ -538,81 +1195,146 @@
           ENDIF
 
-          READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,            &
-                                           highvg_vert, highwsubs_vert,        &
-                                           high_td_lsa_lpt, high_td_lsa_q,     &
-                                           high_td_sub_lpt, high_td_sub_q
-      
-          IF ( ierrn /= 0 )  THEN
-             message_string = 'errors in file LSF_DATA'
-             CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
-          ENDIF
-
-
-          DO  k = nzb, nzt+1
-             IF ( highheight < zu(k) )  THEN
-                lowheight      = highheight
-                lowug_vert     = highug_vert
-                lowvg_vert     = highvg_vert
-                lowwsubs_vert  = highwsubs_vert
-                low_td_lsa_lpt = high_td_lsa_lpt
-                low_td_lsa_q   = high_td_lsa_q
-                low_td_sub_lpt = high_td_sub_lpt
-                low_td_sub_q   = high_td_sub_q
-
-                ierrn = 0
-                READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,      &
-                                                 highvg_vert, highwsubs_vert,  &
-                                                 high_td_lsa_lpt,              &
-                                                 high_td_lsa_q,                &
-                                                 high_td_sub_lpt, high_td_sub_q
-
-                IF ( ierrn /= 0 )  THEN
-                   WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ', &
-                        'is higher than the maximum height in LSF_DATA which ',&
-                        'is ', lowheight, 'm. Interpolation on PALM ',         &
-                        'grid is not possible.'
-                   CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
+!
+!--       Surface values are read in
+          nt     = 0
+          ierrn = 0
+
+          DO WHILE ( time_surf(nt) < end_time )
+             nt = nt + 1
+             READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt),   &
+                                                qsws_surf(nt), pt_surf(nt),    &
+                                                q_surf(nt), p_surf(nt)
+
+             IF ( ierrn /= 0 )  THEN
+               WRITE ( message_string, * ) 'No time dependent surface variables ' //&
+                                 'in&LSF_DATA for end of run found'
+
+                CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
+             ENDIF
+          ENDDO
+
+          IF ( time_surf(1) > end_time )  THEN
+             WRITE ( message_string, * ) 'Time dependent surface variables in ' // &
+                                         '&LSF_DATA set in after end of ' ,        &
+                                         'simulation - lsf_surf is set to FALSE'
+             CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
+             lsf_surf = .FALSE.
+          ENDIF
+
+!
+!--       Go to the end of the list with surface variables
+          DO WHILE ( ierrn == 0 )
+             READ ( finput, *, IOSTAT = ierrn ) r_dummy
+          ENDDO
+
+!
+!--       Profiles of ug, vg and w_subs are read in (large scale forcing)
+
+          nt = 0
+          DO WHILE ( time_vert(nt) < end_time )
+             nt = nt + 1
+             hash = "#"
+             ierrn = 1 ! not zero
+!
+!--          Search for the next line consisting of "# time", 
+!--          from there onwards the profiles will be read
+             DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
+                READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
+                IF ( ierrn < 0 )  THEN 
+                   WRITE( message_string, * ) 'No time dependent vertical profiles',&
+                                    ' in&LSF_DATA for end of run found'
+                   CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
                 ENDIF
-
+             ENDDO
+
+             IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
+
+             READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,&
+                                              lowwsubs_vert, low_td_lsa_lpt,    &
+                                              low_td_lsa_q, low_td_sub_lpt,     &
+                                              low_td_sub_q
+             IF ( ierrn /= 0 )  THEN
+                message_string = 'errors in file LSF_DATA'
+                CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
              ENDIF
 
-!
-!--          Interpolation of prescribed profiles in space 
-             fac = (highheight-zu(k))/(highheight - lowheight)
-
-             ug_vert(k,nt)    = fac * lowug_vert                               &
-                                + ( 1.0_wp - fac ) * highug_vert
-             vg_vert(k,nt)    = fac * lowvg_vert                               &
-                                + ( 1.0_wp - fac ) * highvg_vert
-             wsubs_vert(k,nt) = fac * lowwsubs_vert                            &
-                                + ( 1.0_wp - fac ) * highwsubs_vert
-
-             td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                           &
-                                + ( 1.0_wp - fac ) * high_td_lsa_lpt
-             td_lsa_q(k,nt)   = fac * low_td_lsa_q                             &
-                                + ( 1.0_wp - fac ) * high_td_lsa_q
-             td_sub_lpt(k,nt) = fac * low_td_sub_lpt                           &
-                                + ( 1.0_wp - fac ) * high_td_sub_lpt
-             td_sub_q(k,nt)   = fac * low_td_sub_q                             &
-                                + ( 1.0_wp - fac ) * high_td_sub_q
-
-          ENDDO
-
-       ENDDO 
-
-!
-!--    Large scale vertical velocity has to be zero at the surface
-       wsubs_vert(nzb,:) = 0.0_wp
- 
-       IF ( time_vert(1) > end_time )  THEN
-          WRITE ( message_string, * ) 'Time dependent large scale profile ',   &
-                             'forcing from&LSF_DATA sets in after end of ' ,   &
-                             'simulation - lsf_vert is set to FALSE'
-          CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
-          lsf_vert = .FALSE.
-       ENDIF
-
-       CLOSE( finput )
-
+             READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,         &
+                                              highvg_vert, highwsubs_vert,     &
+                                              high_td_lsa_lpt, high_td_lsa_q,  &
+                                              high_td_sub_lpt, high_td_sub_q
+         
+             IF ( ierrn /= 0 )  THEN
+                message_string = 'errors in file LSF_DATA'
+                CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
+             ENDIF
+
+
+             DO  k = nzb, nzt+1
+                IF ( highheight < zu(k) )  THEN
+                   lowheight      = highheight
+                   lowug_vert     = highug_vert
+                   lowvg_vert     = highvg_vert
+                   lowwsubs_vert  = highwsubs_vert
+                   low_td_lsa_lpt = high_td_lsa_lpt
+                   low_td_lsa_q   = high_td_lsa_q
+                   low_td_sub_lpt = high_td_sub_lpt
+                   low_td_sub_q   = high_td_sub_q
+
+                   ierrn = 0
+                   READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert,    &
+                                                    highvg_vert, highwsubs_vert,&
+                                                    high_td_lsa_lpt,            &
+                                                    high_td_lsa_q,              &
+                                                    high_td_sub_lpt, high_td_sub_q
+
+                   IF ( ierrn /= 0 )  THEN
+                      WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ', &
+                           'is higher than the maximum height in LSF_DATA which ',&
+                           'is ', lowheight, 'm. Interpolation on PALM ',         &
+                           'grid is not possible.'
+                      CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
+                   ENDIF
+
+                ENDIF
+
+!
+!--             Interpolation of prescribed profiles in space 
+                fac = (highheight-zu(k))/(highheight - lowheight)
+
+                ug_vert(k,nt)    = fac * lowug_vert                            &
+                                   + ( 1.0_wp - fac ) * highug_vert
+                vg_vert(k,nt)    = fac * lowvg_vert                            &
+                                   + ( 1.0_wp - fac ) * highvg_vert
+                wsubs_vert(k,nt) = fac * lowwsubs_vert                         &
+                                   + ( 1.0_wp - fac ) * highwsubs_vert
+
+                td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt                        &
+                                   + ( 1.0_wp - fac ) * high_td_lsa_lpt
+                td_lsa_q(k,nt)   = fac * low_td_lsa_q                          &
+                                   + ( 1.0_wp - fac ) * high_td_lsa_q
+                td_sub_lpt(k,nt) = fac * low_td_sub_lpt                        &
+                                   + ( 1.0_wp - fac ) * high_td_sub_lpt
+                td_sub_q(k,nt)   = fac * low_td_sub_q                          &
+                                   + ( 1.0_wp - fac ) * high_td_sub_q
+
+             ENDDO
+
+          ENDDO 
+
+!
+!--       Large scale vertical velocity has to be zero at the surface
+          wsubs_vert(nzb,:) = 0.0_wp
+    
+          IF ( time_vert(1) > end_time )  THEN
+             WRITE ( message_string, * ) 'Time dependent large scale profile ',&
+                                'forcing from&LSF_DATA sets in after end of ' ,&
+                                'simulation - lsf_vert is set to FALSE'
+             CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
+             lsf_vert = .FALSE.
+          ENDIF
+
+          CLOSE( finput )
+
+       ENDIF
 
     END SUBROUTINE lsf_init
@@ -1396,3 +2118,4 @@
     END SUBROUTINE nudge_ref
 
+
  END MODULE lsf_nudging_mod
