Index: /palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 146)
+++ /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 147)
@@ -15,7 +15,13 @@
 necessary for the current subdomain.
 
+The turbulence recycling method of Kataoka and Mizuno (2002) has been
+implemented. A pre-run with smaller model domain can be used to initialize
+the real run, in order to spare the CPU time needed to get the turbulence
+into a quasi-stationary state (see initializing_actions =
+'read_data_for_recycling').
+
 User-defined spectra.
 
-calc_spectra, check_open, data_output_spectra, init_pegrid, modules, netcdf, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary
+calc_spectra, check_open, check_parameters, data_output_spectra, init_pegrid, init_3d_model, modules, netcdf, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary
 
 
@@ -50,5 +56,7 @@
 contained uneccessary time levels. (read_3d_binary, write_3d_binary)
 Bugfix: extra '*' removed in user_statistics sample code (user_interface)
+Bugfix: a stop command was missing in some cases of the parallel branch (local_stop)
 
-flow_statistics, plant_canopy_model, read_3d_binary, user_interface, write_3d_binary
 
+flow_statistics, local_stop, plant_canopy_model, read_3d_binary, user_interface, write_3d_binary
+
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 146)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 147)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Case of reading data for recycling included in initializing_actions 
 !
 ! Former revisions:
@@ -506,5 +506,6 @@
     ENDIF
 
-    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
+         TRIM( initializing_actions ) /= 'read_data_for_recycling' )  THEN
 !
 !--    No model continuation run; several initialising actions are possible
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 146)
+++ /palm/trunk/SOURCE/header.f90	(revision 147)
@@ -108,12 +108,12 @@
     IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
        run_classification = '3D - restart run'
+    ELSEIF ( TRIM( initializing_actions ) == 'read_data_for_recycling' )  THEN
+       run_classification = '3D - run using 3D - prerun data'
+    ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 )  THEN
+       run_classification = '3D - run without 1D - prerun'
+    ELSEIF ( INDEX(initializing_actions, 'set_1d-model_profiles') /= 0 ) THEN
+       run_classification = '3D - run with 1D - prerun'
     ELSE
-       IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 )  THEN
-          run_classification = '3D - run without 1D - prerun'
-       ELSEIF ( INDEX(initializing_actions, 'set_1d-model_profiles') /= 0 ) THEN
-          run_classification = '3D - run with 1D - prerun'
-       ELSE
-          PRINT*,'+++ header:  unknown action(s): ',initializing_actions
-       ENDIF
+       PRINT*,'+++ header:  unknown action(s): ',initializing_actions
     ENDIF
     IF ( ocean )  THEN
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 146)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 147)
@@ -8,5 +8,6 @@
 ! -----------------
 ! Allocation of hom_sum moved to parin, initialization of spectrum_x|y directly
-! after allocating theses arrays
+! after allocating theses arrays,
+! read data for recycling added as new initialization option
 !
 ! Former revisions:
@@ -366,5 +367,6 @@
 !
 !-- Initialize model variables
-    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
+         TRIM( initializing_actions ) /= 'read_data_for_recycling' )  THEN
 !
 !--    First model run of a possible job queue.
@@ -907,9 +909,26 @@
        ENDIF
 
-    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' ) &
+    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.    &
+             TRIM( initializing_actions ) == 'read_data_for_recycling' )  &
     THEN
 !
+!--    When reading data for initializing the recycling method, first read
+!--    some of the global variables from restart file
+       IF ( TRIM( initializing_actions ) == 'read_data_for_recycling' )  THEN
+          WRITE (9,*) 'before read_parts_of_var_list'
+          CALL local_flush( 9 )
+          CALL read_parts_of_var_list
+          WRITE (9,*) 'after read_parts_of_var_list'
+          CALL local_flush( 9 )
+          CALL close_file( 13 )
+       ENDIF
+
+!
 !--    Read binary data from restart file
+          WRITE (9,*) 'before read_3d_binary'
+          CALL local_flush( 9 )
        CALL read_3d_binary
+          WRITE (9,*) 'after read_3d_binary'
+          CALL local_flush( 9 )
 
 !
Index: /palm/trunk/SOURCE/local_stop.f90
===================================================================
--- /palm/trunk/SOURCE/local_stop.f90	(revision 146)
+++ /palm/trunk/SOURCE/local_stop.f90	(revision 147)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Bugfix: a stop command was missing in some cases of the parallel branch
 ! 
 !
@@ -34,4 +34,5 @@
     IF ( coupling_mode == 'uncoupled' )  THEN
        CALL MPI_FINALIZE( ierr )
+       STOP
     ELSE
 
@@ -56,4 +57,5 @@
              ENDIF
              CALL MPI_FINALIZE( ierr )
+             STOP
 
           CASE ( 1 )
@@ -64,4 +66,5 @@
              ENDIF
              CALL MPI_FINALIZE( ierr )
+             STOP
 
           CASE ( 2 )
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 146)
+++ /palm/trunk/SOURCE/modules.f90	(revision 147)
@@ -5,5 +5,6 @@
 ! Actual revisions:
 ! -----------------
-! +hor_index_bounds, hor_index_bounds_previous_run, numprocs_previous_run
+! +hor_index_bounds, hor_index_bounds_previous_run, numprocs_previous_run,
+! nx_on_file, ny_on_file
 ! -myid_char_14
 !
@@ -595,7 +596,7 @@
 !------------------------------------------------------------------------------!
 
-    INTEGER ::  ngp_sums, nnx, nx = 0, nxa, nxl, nxlu, nxr, nxra, nny, ny = 0, &
-                nya, nyn, nyna, nys, nysv, nnz, nz = 0, nza, nzb, nzb_diff,    &
-                nzt, nzta, nzt_diff
+    INTEGER ::  ngp_sums, nnx, nx = 0, nxa, nxl, nxlu, nxr, nxra, nx_on_file,  &
+                nny, ny = 0, nya, nyn, nyna, nys, nysv, ny_on_file, nnz,       &
+                nz = 0, nza, nzb, nzb_diff, nzt, nzta, nzt_diff
 
     INTEGER, DIMENSION(:), ALLOCATABLE ::                                      &
Index: /palm/trunk/SOURCE/parin.f90
===================================================================
--- /palm/trunk/SOURCE/parin.f90	(revision 146)
+++ /palm/trunk/SOURCE/parin.f90	(revision 147)
@@ -231,5 +231,7 @@
           CALL local_stop
        ENDIF
-
+!
+!--    ATTENTION: in case of changes to the following statement please also
+!--    check the allocate statement in routine read_var_list
        ALLOCATE( lad(0:nz+1),pt_init(0:nz+1), q_init(0:nz+1), sa_init(0:nz+1), &
                  ug(0:nz+1), u_init(0:nz+1), v_init(0:nz+1), vg(0:nz+1),       &
Index: /palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/read_3d_binary.f90	(revision 146)
+++ /palm/trunk/SOURCE/read_3d_binary.f90	(revision 147)
@@ -5,5 +5,7 @@
 ! -----------------
 ! Files from which restart data are to be read are determined and subsequently
-! opened,
+! opened. The total domain on the restart file is allowed to be smaller than
+! the current total domain. In this case it will be periodically mapped on the
+! current domain (needed for recycling method).
 ! +call of user_read_restart_data, -dopr_time_count,
 ! hom_sum, volume_flow_area, volume_flow_initial moved to read_var_list,
@@ -62,8 +64,11 @@
 
     INTEGER ::  files_to_be_opened, i, idum1, j, myid_on_file,                &
-                numprocs_on_file, nxlc, nxl_on_file, nxrc, nxr_on_file, nync, &
-                nyn_on_file, nysc, nys_on_file, nzb_on_file, nzt_on_file
-
-    INTEGER, DIMENSION(:), ALLOCATABLE ::  file_list
+                numprocs_on_file, nxlc, nxlf, nxlpr, nxl_on_file, nxrc, nxrf, &
+                nxrpr, nxr_on_file, nync, nynf, nynpr, nyn_on_file, nysc,     &
+                nysf, nyspr, nys_on_file, nzb_on_file, nzt_on_file, offset_x, &
+                offset_y
+
+    INTEGER, DIMENSION(numprocs_previous_run*4) ::  file_list, nxlfa, nxrfa, &
+                nynfa, nysfa, offset_xa, offset_ya
 
     REAL, DIMENSION(:,:), ALLOCATABLE     ::  tmp_2d
@@ -80,24 +85,187 @@
 !-- of this PE
     files_to_be_opened = 0
-    ALLOCATE( file_list(numprocs_previous_run) )
 
     DO  i = 1, numprocs_previous_run
 
+       nxlpr = hor_index_bounds_previous_run(1,i-1)
+       nxrpr = hor_index_bounds_previous_run(2,i-1)
+       nyspr = hor_index_bounds_previous_run(3,i-1)
+       nynpr = hor_index_bounds_previous_run(4,i-1)
+
+!
+!--    Determine the offsets. They may be non-zero in case that the total domain
+!--    on file is smaller than the current total domain.
+       offset_x = ( nxl / ( nx_on_file + 1 ) ) * ( nx_on_file + 1 )
+       offset_y = ( nys / ( ny_on_file + 1 ) ) * ( ny_on_file + 1 )
+
 !
 !--    Only data which overlap with the current subdomain have to be read
-       IF ( hor_index_bounds_previous_run(1,i-1) <= nxr  .AND.  &
-            hor_index_bounds_previous_run(2,i-1) >= nxl  .AND.  &
-            hor_index_bounds_previous_run(3,i-1) <= nyn  .AND.  &
-            hor_index_bounds_previous_run(4,i-1) >= nys )  THEN
+       IF ( nxlpr <= nxr-offset_x  .AND.  nxrpr >= nxl-offset_x  .AND.  &
+            nyspr <= nyn-offset_y  .AND.  nynpr >= nys-offset_y )  THEN
 
           files_to_be_opened = files_to_be_opened + 1
           file_list(files_to_be_opened) = i-1
+!
+!--       Index bounds of overlapping data
+          nxlfa(files_to_be_opened) = MAX( nxl-offset_x, nxlpr )
+          nxrfa(files_to_be_opened) = MIN( nxr-offset_x, nxrpr )
+          nysfa(files_to_be_opened) = MAX( nys-offset_y, nyspr )
+          nynfa(files_to_be_opened) = MIN( nyn-offset_y, nynpr )
+
+          WRITE (9,*) '*** reading from file: ', i
+          WRITE (9,*) '    index bounds on file:'
+          WRITE (9,*) '       nxlpr=', nxlpr, ' nxrpr=', nxrpr
+          WRITE (9,*) '       nyspr=', nyspr, ' nynpr=', nynpr
+          WRITE (9,*) '    index bounds of current subdmain:'
+          WRITE (9,*) '       nxl  =', nxl, ' nxr  =', nxr
+          WRITE (9,*) '       nys  =', nys, ' nyn  =', nyn
+          WRITE (9,*) '    offset used:'
+          WRITE (9,*) '       offset_x=', offset_x, ' offset_y=', offset_y
+          WRITE (9,*) '    bounds of overlapping data:'
+          WRITE (9,*) '       nxlfa=', nxlfa(files_to_be_opened), ' nxrfa=', nxrfa(files_to_be_opened)
+          WRITE (9,*) '       nysfa=', nysfa(files_to_be_opened), ' nynfa=', nynfa(files_to_be_opened)
+          CALL local_flush( 9 )
+          offset_xa(files_to_be_opened) = offset_x
+          offset_ya(files_to_be_opened) = offset_y
 
        ENDIF
 
+!
+!--    If the total domain on file is smaller than the current total domain,
+!--    and if the current subdomain extends beyond the limits of the total
+!--    domain of file, the respective file may be opened again (three times
+!--    maximum) to read the still missing parts, which are then added
+!--    "cyclically".
+!--    Overlap along x:
+       IF ( ( nxr - offset_x ) > nx_on_file )  THEN
+
+          offset_x = offset_x + ( nx_on_file + 1 )
+
+          IF ( nxlpr <= nxr-offset_x  .AND.  nxrpr >= nxl-offset_x  .AND.  &
+               nyspr <= nyn-offset_y  .AND.  nynpr >= nys-offset_y )  THEN
+
+             files_to_be_opened = files_to_be_opened + 1
+             file_list(files_to_be_opened) = i-1
+!
+!--          Index bounds of overlapping data
+             nxlfa(files_to_be_opened) = MAX( nxl-offset_x, nxlpr )
+             nxrfa(files_to_be_opened) = MIN( nxr-offset_x, nxrpr )
+             nysfa(files_to_be_opened) = MAX( nys-offset_y, nyspr )
+             nynfa(files_to_be_opened) = MIN( nyn-offset_y, nynpr )
+
+          WRITE (9,*) '*** reading from file: ', i
+          WRITE (9,*) '    index bounds on file:'
+          WRITE (9,*) '       nxlpr=', nxlpr, ' nxrpr=', nxrpr
+          WRITE (9,*) '       nyspr=', nyspr, ' nynpr=', nynpr
+          WRITE (9,*) '    index bounds of current subdmain:'
+          WRITE (9,*) '       nxl  =', nxl, ' nxr  =', nxr
+          WRITE (9,*) '       nys  =', nys, ' nyn  =', nyn
+          WRITE (9,*) '    offset used:'
+          WRITE (9,*) '       offset_x=', offset_x, ' offset_y=', offset_y
+          WRITE (9,*) '    bounds of overlapping data:'
+          WRITE (9,*) '       nxlfa=', nxlfa(files_to_be_opened), ' nxrfa=', nxrfa(files_to_be_opened)
+          WRITE (9,*) '       nysfa=', nysfa(files_to_be_opened), ' nynfa=', nynfa(files_to_be_opened)
+          CALL local_flush( 9 )
+             offset_xa(files_to_be_opened) = offset_x
+             offset_ya(files_to_be_opened) = offset_y
+
+          ENDIF
+
+          offset_x = offset_x - ( nx_on_file + 1 )
+
+       ENDIF
+
+
+!
+!--    Overlap along y:
+       IF ( ( nyn - offset_y ) > ny_on_file )  THEN
+
+          offset_y = offset_y + ( ny_on_file + 1 )
+
+          IF ( nxlpr <= nxr-offset_x  .AND.  nxrpr >= nxl-offset_x  .AND.  &
+               nyspr <= nyn-offset_y  .AND.  nynpr >= nys-offset_y )  THEN
+
+             files_to_be_opened = files_to_be_opened + 1
+             file_list(files_to_be_opened) = i-1
+!
+!--          Index bounds of overlapping data
+             nxlfa(files_to_be_opened) = MAX( nxl-offset_x, nxlpr )
+             nxrfa(files_to_be_opened) = MIN( nxr-offset_x, nxrpr )
+             nysfa(files_to_be_opened) = MAX( nys-offset_y, nyspr )
+             nynfa(files_to_be_opened) = MIN( nyn-offset_y, nynpr )
+
+          WRITE (9,*) '*** reading from file: ', i
+          WRITE (9,*) '    index bounds on file:'
+          WRITE (9,*) '       nxlpr=', nxlpr, ' nxrpr=', nxrpr
+          WRITE (9,*) '       nyspr=', nyspr, ' nynpr=', nynpr
+          WRITE (9,*) '    index bounds of current subdmain:'
+          WRITE (9,*) '       nxl  =', nxl, ' nxr  =', nxr
+          WRITE (9,*) '       nys  =', nys, ' nyn  =', nyn
+          WRITE (9,*) '    offset used:'
+          WRITE (9,*) '       offset_x=', offset_x, ' offset_y=', offset_y
+          WRITE (9,*) '    bounds of overlapping data:'
+          WRITE (9,*) '       nxlfa=', nxlfa(files_to_be_opened), ' nxrfa=', nxrfa(files_to_be_opened)
+          WRITE (9,*) '       nysfa=', nysfa(files_to_be_opened), ' nynfa=', nynfa(files_to_be_opened)
+          CALL local_flush( 9 )
+             offset_xa(files_to_be_opened) = offset_x
+             offset_ya(files_to_be_opened) = offset_y
+
+          ENDIF
+
+          offset_y = offset_y - ( ny_on_file + 1 )
+
+       ENDIF
+
+!--    Overlap along x and y:
+       IF ( ( nxr - offset_x ) > nx_on_file  .AND.  &
+            ( nyn - offset_y ) > ny_on_file )  THEN
+
+          offset_x = offset_x + ( nx_on_file + 1 )
+          offset_y = offset_y + ( ny_on_file + 1 )
+
+          IF ( nxlpr <= nxr-offset_x  .AND.  nxrpr >= nxl-offset_x  .AND.  &
+               nyspr <= nyn-offset_y  .AND.  nynpr >= nys-offset_y )  THEN
+
+             files_to_be_opened = files_to_be_opened + 1
+             file_list(files_to_be_opened) = i-1
+!
+!--          Index bounds of overlapping data
+             nxlfa(files_to_be_opened) = MAX( nxl-offset_x, nxlpr )
+             nxrfa(files_to_be_opened) = MIN( nxr-offset_x, nxrpr )
+             nysfa(files_to_be_opened) = MAX( nys-offset_y, nyspr )
+             nynfa(files_to_be_opened) = MIN( nyn-offset_y, nynpr )
+
+          WRITE (9,*) '*** reading from file: ', i
+          WRITE (9,*) '    index bounds on file:'
+          WRITE (9,*) '       nxlpr=', nxlpr, ' nxrpr=', nxrpr
+          WRITE (9,*) '       nyspr=', nyspr, ' nynpr=', nynpr
+          WRITE (9,*) '    index bounds of current subdmain:'
+          WRITE (9,*) '       nxl  =', nxl, ' nxr  =', nxr
+          WRITE (9,*) '       nys  =', nys, ' nyn  =', nyn
+          WRITE (9,*) '    offset used:'
+          WRITE (9,*) '       offset_x=', offset_x, ' offset_y=', offset_y
+          WRITE (9,*) '    bounds of overlapping data:'
+          WRITE (9,*) '       nxlfa=', nxlfa(files_to_be_opened), ' nxrfa=', nxrfa(files_to_be_opened)
+          WRITE (9,*) '       nysfa=', nysfa(files_to_be_opened), ' nynfa=', nynfa(files_to_be_opened)
+          CALL local_flush( 9 )
+             offset_xa(files_to_be_opened) = offset_x
+             offset_ya(files_to_be_opened) = offset_y
+
+          ENDIF
+
+          offset_x = offset_x - ( nx_on_file + 1 )
+          offset_y = offset_y - ( ny_on_file + 1 )
+
+       ENDIF
+
     ENDDO
 
+!
+!-- Save the id-string of the current process, since myid_char may now be used
+!-- to open files created by PEs with other id.
     myid_char_save = myid_char
 
+!
+!-- Test output (remove later)
     DO i = 1, numprocs_previous_run
        WRITE (9,*) 'i=',i-1, ' ibs= ',hor_index_bounds_previous_run(1:4,i-1)
@@ -128,5 +296,9 @@
 !--    first.
        CALL check_open( 13 )
+       WRITE (9,*) 'before skipping'
+       CALL local_flush( 9 )
        IF ( j == 0 )  CALL skip_var_list
+       WRITE (9,*) 'skipping done'
+       CALL local_flush( 9 )
 
 !
@@ -228,10 +400,14 @@
 
 !
-!--    Determine the index range of those gridpoints to be copied from
-!--    the subdomains on the restart files to the current subdomain
-       nxlc = MAX( nxl, nxl_on_file )
-       nxrc = MIN( nxr, nxr_on_file )
-       nysc = MAX( nys, nys_on_file )
-       nync = MIN( nyn, nyn_on_file )
+!--    Get the index range of the subdomain on file which overlap with the
+!--    current subdomain
+       nxlf = nxlfa(i)
+       nxlc = nxlfa(i) + offset_xa(i)
+       nxrf = nxrfa(i)
+       nxrc = nxrfa(i) + offset_xa(i)
+       nysf = nysfa(i)
+       nysc = nysfa(i) + offset_ya(i)
+       nynf = nynfa(i)
+       nync = nynfa(i) + offset_ya(i)
 
 !
@@ -244,4 +420,6 @@
        DO  WHILE ( TRIM( field_chr ) /= '*** end ***' )
 
+          WRITE (9,*) 'var = ', field_chr
+          CALL local_flush( 9 )
           SELECT CASE ( TRIM( field_chr ) )
 
@@ -249,5 +427,5 @@
                 READ ( 13 )  tmp_3d
                 e(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'e_av' )
@@ -257,10 +435,10 @@
                 READ ( 13 )  tmp_3d
                 e_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'e_m' )
                 READ ( 13 )  tmp_3d
                 e_m(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'iran' ) ! matching random numbers is still unresolved issue
@@ -270,20 +448,20 @@
                 READ ( 13 )  tmp_3d
                 kh(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'kh_m' )
                 READ ( 13 )  tmp_3d
                 kh_m(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'km' )
                 READ ( 13 )  tmp_3d
                 km(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'km_m' )
                 READ ( 13 )  tmp_3d
                 km_m(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'lwp_av' )
@@ -293,10 +471,10 @@
                 READ ( 13 )  tmp_2d
                 lwp_av(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'p' )
                 READ ( 13 )  tmp_3d
                 p(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'p_av' )
@@ -306,5 +484,5 @@
                 READ ( 13 )  tmp_3d
                 p_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'pc_av' )
@@ -314,5 +492,5 @@
                 READ ( 13 )  tmp_3d
                 pc_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'pr_av' )
@@ -322,10 +500,10 @@
                 READ ( 13 )  tmp_3d
                 pr_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'precipitation_amount' )
                 READ ( 13 )  tmp_2d
                 precipitation_amount(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'precipitation_rate_a' )
@@ -335,10 +513,10 @@
                 READ ( 13 )  tmp_2d
                 precipitation_rate_av(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'pt' )
                 READ ( 13 )  tmp_3d
                 pt(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'pt_av' )
@@ -348,15 +526,15 @@
                 READ ( 13 )  tmp_3d
                 pt_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'pt_m' )
                 READ ( 13 )  tmp_3d
                 pt_m(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'q' )
                 READ ( 13 )  tmp_3d
                 q(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'q_av' )
@@ -366,15 +544,15 @@
                 READ ( 13 )  tmp_3d
                 q_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'q_m' )
                 READ ( 13 )  tmp_3d
                 q_m(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'ql' )
                 READ ( 13 )  tmp_3d
                 ql(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'ql_av' )
@@ -384,5 +562,5 @@
                 READ ( 13 )  tmp_3d
                 ql_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'ql_c_av' )
@@ -392,5 +570,5 @@
                 READ ( 13 )  tmp_3d
                 ql_c_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'ql_v_av' )
@@ -400,5 +578,5 @@
                 READ ( 13 )  tmp_3d
                 ql_v_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'ql_vp_av' )
@@ -408,30 +586,30 @@
                 READ ( 13 )  tmp_3d
                 ql_vp_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'qs' )
                 READ ( 13 )  tmp_2d
                 qs(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'qsws' )
                 READ ( 13 )  tmp_2d
                 qsws(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'qsws_m' )
                 READ ( 13 )  tmp_2d
                 qsws_m(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'qswst' )
                 READ ( 13 )  tmp_2d
                 qswst(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'qswst_m' )
                 READ ( 13 )  tmp_2d
                 qswst_m(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'qv_av' )
@@ -441,5 +619,5 @@
                 READ ( 13 )  tmp_3d
                 qv_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'random_iv' )  ! still unresolved issue
@@ -453,15 +631,15 @@
                 READ ( 13 )  tmp_3d
                 rho_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'rif' )
                 READ ( 13 )  tmp_2d
                 rif(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'rif_m' )
                 READ ( 13 )  tmp_2d
                 rif_m(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'rif_wall' )
@@ -470,5 +648,5 @@
                 READ ( 13 )  tmp_4d
                 rif_wall(:,nysc-1:nync+1,nxlc-1:nxrc+1,:) = &
-                                         tmp_4d(:,nysc-1:nync+1,nxlc-1:nxrc+1,:)
+                                         tmp_4d(:,nysf-1:nynf+1,nxlf-1:nxrf+1,:)
                 DEALLOCATE( tmp_4d )
 
@@ -479,10 +657,10 @@
                 READ ( 13 )  tmp_3d
                 s_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'sa' )
                 READ ( 13 )  tmp_3d
                 sa(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'sa_av' )
@@ -492,25 +670,25 @@
                 READ ( 13 )  tmp_3d
                 sa_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'saswsb' )
                 READ ( 13 )  tmp_2d
                 saswsb(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'saswst' )
                 READ ( 13 )  tmp_2d
                 saswst(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'shf' )
                 READ ( 13 )  tmp_2d
                 shf(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'shf_m' )
                 READ ( 13 )  tmp_2d
                 shf_m(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'spectrum_x' )
@@ -523,5 +701,5 @@
                 READ ( 13 )  tmp_2d
                 ts(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'ts_av' )
@@ -531,20 +709,20 @@
                 READ ( 13 )  tmp_2d
                 ts_av(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'tswst' )
                 READ ( 13 )  tmp_2d
                 tswst(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'tswst_m' )
                 READ ( 13 )  tmp_2d
                 tswst_m(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'u' )
                 READ ( 13 )  tmp_3d
                 u(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'u_av' )
@@ -554,10 +732,10 @@
                 READ ( 13 )  tmp_3d
                 u_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'u_m' )
                 READ ( 13 )  tmp_3d
                 u_m(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'u_m_l' )
@@ -565,5 +743,5 @@
                 READ ( 13 )  tmp_3dw
                 IF ( outflow_l )  THEN
-                   u_m_l(:,nysc-1:nync+1,:) = tmp_3dw(:,nysc-1:nync+1,:)
+                   u_m_l(:,nysc-1:nync+1,:) = tmp_3dw(:,nysf-1:nynf+1,:)
                 ENDIF
                 DEALLOCATE( tmp_3dw )
@@ -574,5 +752,5 @@
                 READ ( 13 )  tmp_3dw
                 IF ( outflow_n )  THEN
-                   u_m_n(:,:,nxlc-1:nxrc+1) = tmp_3dw(:,:,nxlc-1:nxrc+1)
+                   u_m_n(:,:,nxlc-1:nxrc+1) = tmp_3dw(:,:,nxlf-1:nxrf+1)
                 ENDIF
                 DEALLOCATE( tmp_3dw )
@@ -583,5 +761,5 @@
                 READ ( 13 )  tmp_3dw
                 IF ( outflow_r )  THEN
-                   u_m_r(:,nysc-1:nync+1,:) = tmp_3dw(:,nysc-1:nync+1,:)
+                   u_m_r(:,nysc-1:nync+1,:) = tmp_3dw(:,nysf-1:nynf+1,:)
                 ENDIF
                 DEALLOCATE( tmp_3dw )
@@ -592,5 +770,5 @@
                 READ ( 13 )  tmp_3dw
                 IF ( outflow_s )  THEN
-                   u_m_s(:,:,nxlc-1:nxrc+1) = tmp_3dw(:,:,nxlc-1:nxrc+1)
+                   u_m_s(:,:,nxlc-1:nxrc+1) = tmp_3dw(:,:,nxlf-1:nxrf+1)
                 ENDIF
                 DEALLOCATE( tmp_3dw )
@@ -599,25 +777,25 @@
                 READ ( 13 )  tmp_2d
                 us(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'usws' )
                 READ ( 13 )  tmp_2d
                 usws(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'uswst' )
                 READ ( 13 )  tmp_2d
                 uswst(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'usws_m' )
                 READ ( 13 )  tmp_2d
                 usws_m(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'uswst_m' )
                 READ ( 13 )  tmp_2d
                 uswst_m(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'us_av' )
@@ -627,10 +805,10 @@
                 READ ( 13 )  tmp_2d
                 us_av(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'v' )
                 READ ( 13 )  tmp_3d
                 v(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'v_av' )
@@ -640,10 +818,10 @@
                 READ ( 13 )  tmp_3d
                 v_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'v_m' )
                 READ ( 13 )  tmp_3d
                 v_m(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'v_m_l' )
@@ -651,5 +829,5 @@
                 READ ( 13 )  tmp_3dw
                 IF ( outflow_l )  THEN
-                   v_m_l(:,nysc-1:nync+1,:) = tmp_3dw(:,nysc-1:nync+1,:)
+                   v_m_l(:,nysc-1:nync+1,:) = tmp_3dw(:,nysf-1:nynf+1,:)
                 ENDIF
                 DEALLOCATE( tmp_3dw )
@@ -660,5 +838,5 @@
                 READ ( 13 )  tmp_3dw
                 IF ( outflow_n )  THEN
-                   v_m_n(:,:,nxlc-1:nxrc+1) = tmp_3dw(:,:,nxlc-1:nxrc+1)
+                   v_m_n(:,:,nxlc-1:nxrc+1) = tmp_3dw(:,:,nxlf-1:nxrf+1)
                 ENDIF
                 DEALLOCATE( tmp_3dw )
@@ -669,5 +847,5 @@
                 READ ( 13 )  tmp_3dw
                 IF ( outflow_r )  THEN
-                   v_m_r(:,nysc-1:nync+1,:) = tmp_3dw(:,nysc-1:nync+1,:)
+                   v_m_r(:,nysc-1:nync+1,:) = tmp_3dw(:,nysf-1:nynf+1,:)
                 ENDIF
                 DEALLOCATE( tmp_3dw )
@@ -678,5 +856,5 @@
                 READ ( 13 )  tmp_3dw
                 IF ( outflow_s )  THEN
-                   v_m_s(:,:,nxlc-1:nxrc+1) = tmp_3dw(:,:,nxlc-1:nxrc+1)
+                   v_m_s(:,:,nxlc-1:nxrc+1) = tmp_3dw(:,:,nxlf-1:nxrf+1)
                 ENDIF
                 DEALLOCATE( tmp_3dw )
@@ -685,5 +863,5 @@
                 READ ( 13 )  tmp_3d
                 vpt(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'vpt_av' )
@@ -693,35 +871,35 @@
                 READ ( 13 )  tmp_3d
                 vpt_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'vpt_m' )
                 READ ( 13 )  tmp_3d
                 vpt_m(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'vsws' )
                 READ ( 13 )  tmp_2d
                 vsws(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'vswst' )
                 READ ( 13 )  tmp_2d
                 vswst(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'vsws_m' )
                 READ ( 13 )  tmp_2d
                 vsws_m(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'vswst_m' )
                 READ ( 13 )  tmp_2d
                 vswst_m(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'w' )
                 READ ( 13 )  tmp_3d
                 w(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'w_av' )
@@ -731,10 +909,10 @@
                 READ ( 13 )  tmp_3d
                 w_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'w_m' )
                 READ ( 13 )  tmp_3d
                 w_m(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_3d(:,nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'w_m_l' )
@@ -742,5 +920,5 @@
                 READ ( 13 )  tmp_3dw
                 IF ( outflow_l )  THEN
-                   w_m_l(:,nysc-1:nync+1,:) = tmp_3dw(:,nysc-1:nync+1,:)
+                   w_m_l(:,nysc-1:nync+1,:) = tmp_3dw(:,nysf-1:nynf+1,:)
                 ENDIF
                 DEALLOCATE( tmp_3dw )
@@ -751,5 +929,5 @@
                 READ ( 13 )  tmp_3dw
                 IF ( outflow_n )  THEN
-                   w_m_n(:,:,nxlc-1:nxrc+1) = tmp_3dw(:,:,nxlc-1:nxrc+1)
+                   w_m_n(:,:,nxlc-1:nxrc+1) = tmp_3dw(:,:,nxlf-1:nxrf+1)
                 ENDIF
                 DEALLOCATE( tmp_3dw )
@@ -760,5 +938,5 @@
                 READ ( 13 )  tmp_3dw
                 IF ( outflow_r )  THEN
-                   w_m_r(:,nysc-1:nync+1,:) = tmp_3dw(:,nysc-1:nync+1,:)
+                   w_m_r(:,nysc-1:nync+1,:) = tmp_3dw(:,nysf-1:nynf+1,:)
                 ENDIF
                 DEALLOCATE( tmp_3dw )
@@ -769,5 +947,5 @@
                 READ ( 13 )  tmp_3dw
                 IF ( outflow_s )  THEN
-                   w_m_s(:,:,nxlc-1:nxrc+1) = tmp_3dw(:,:,nxlc-1:nxrc+1)
+                   w_m_s(:,:,nxlc-1:nxrc+1) = tmp_3dw(:,:,nxlf-1:nxrf+1)
                 ENDIF
                 DEALLOCATE( tmp_3dw )
@@ -776,5 +954,5 @@
                 READ ( 13 )  tmp_2d
                 z0(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE ( 'z0_av' )
@@ -784,5 +962,5 @@
                 READ ( 13 )  tmp_2d
                 z0_av(nysc-1:nync+1,nxlc-1:nxrc+1) = &
-                                          tmp_2d(nysc-1:nync+1,nxlc-1:nxrc+1)
+                                          tmp_2d(nysf-1:nynf+1,nxlf-1:nxrf+1)
 
              CASE DEFAULT
@@ -801,5 +979,7 @@
 !
 !--    Read user-defined restart data
-       CALL user_read_restart_data
+       CALL user_read_restart_data( nxlc, nxlf, nxl_on_file, nxrc, nxrf, &
+                                    nxr_on_file, nync, nynf, nyn_on_file, &
+                                    nysc, nysf, nys_on_file, tmp_2d, tmp_3d )
 
 !
Index: /palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/read_var_list.f90	(revision 146)
+++ /palm/trunk/SOURCE/read_var_list.f90	(revision 147)
@@ -8,5 +8,5 @@
 ! hom_sum, volume_flow_area, volume_flow_initial moved from
 ! read_3d_binary to here,
-! routine skip_var_list added at the end
+! routines read_parts_of_var_list and skip_var_list added at the end
 !
 ! Former revisions:
@@ -148,6 +148,7 @@
     READ ( 13 )  statistic_regions
     IF ( .NOT. ALLOCATED( ug ) )  THEN
-       ALLOCATE( ug(0:nz+1), u_init(0:nz+1), vg(0:nz+1), v_init(0:nz+1), & 
-                 pt_init(0:nz+1), q_init(0:nz+1), sa_init(0:nz+1),       &
+       ALLOCATE( lad(0:nz+1), ug(0:nz+1), u_init(0:nz+1), vg(0:nz+1),    &
+                 v_init(0:nz+1), pt_init(0:nz+1), q_init(0:nz+1),        &
+                 sa_init(0:nz+1),                                        &
                  hom(0:nz+1,2,pr_palm+max_pr_user,0:statistic_regions),  &
                  hom_sum(0:nz+1,pr_palm+max_pr_user,0:statistic_regions) )
@@ -314,6 +315,8 @@
           CASE ( 'nx' )
              READ ( 13 )  nx
+             nx_on_file = nx
           CASE ( 'ny' )
              READ ( 13 )  ny
+             ny_on_file = ny
           CASE ( 'ocean' )
              READ ( 13 )  ocean
@@ -542,4 +545,174 @@
 
 
+ SUBROUTINE read_parts_of_var_list
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! Skipping the global control variables from restart-file (binary format)
+! except some informations needed when reading restart data from a previous
+! run which used a smaller total domain.
+!------------------------------------------------------------------------------!
+
+    USE control_parameters
+    USE indices
+    USE pegrid
+    USE statistics
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=10) ::  version_on_file
+    CHARACTER (LEN=30) ::  variable_chr
+
+    INTEGER ::  idum, max_pr_user_on_file, nz_on_file, &
+                statistic_regions_on_file
+
+    REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  hom_on_file
+
+
+    CALL check_open( 13 )
+
+    WRITE (9,*) 'rpovl: after check open 13'
+    CALL local_flush( 9 )
+    READ ( 13 )  version_on_file
+
+!
+!-- Read number of PEs and horizontal index bounds of all PEs used in previous
+!-- run
+    READ ( 13 )  variable_chr
+    IF ( TRIM( variable_chr ) /= 'numprocs' )  THEN
+       PRINT*, '+++ read_parts_of_var_list: numprocs not found in data from ', &
+                    'prior run on PE ', myid
+       CALL local_stop
+    ENDIF
+    READ ( 13 )  numprocs_previous_run
+
+    IF ( .NOT. ALLOCATED( hor_index_bounds_previous_run ) )  THEN
+       ALLOCATE( hor_index_bounds_previous_run(4,0:numprocs_previous_run-1) )
+    ENDIF
+
+    READ ( 13 )  variable_chr
+    IF ( TRIM( variable_chr ) /= 'hor_index_bounds' )  THEN
+       PRINT*, '+++ read_parts_of_var_list: hor_index_bounds not found in da', &
+                    'ta from prior run on PE ', myid
+       CALL local_stop
+    ENDIF
+    READ ( 13 )  hor_index_bounds_previous_run
+
+!
+!-- Read vertical number of gridpoints and number of different areas used
+!-- for computing statistics. Allocate arrays depending on these values,
+!-- which are needed for the following read instructions.
+    READ ( 13 )  variable_chr
+    IF ( TRIM( variable_chr ) /= 'nz' )  THEN
+       PRINT*, '+++ read_parts_of_var_list: nz not found in restart data file'
+       CALL local_stop
+    ENDIF
+    READ ( 13 )  nz_on_file
+    IF ( nz_on_file /= nz )  THEN
+       IF ( myid == 0 )  THEN
+          PRINT*, '+++ read_parts_of_var_list: mismatch concerning number of', &
+                       ' gridpoints along z'
+          PRINT*, '                   nz on file    = "', nz_on_file, '"'
+          PRINT*, '                   nz from run   = "', nz, '"'
+       ENDIF
+       CALL local_stop
+    ENDIF
+
+    READ ( 13 )  variable_chr
+    IF ( TRIM( variable_chr ) /= 'max_pr_user' )  THEN
+       PRINT*, '+++ read_parts_of_var_list: max_pr_user not found in restart', &
+                    ' data file'
+       CALL local_stop
+    ENDIF
+    READ ( 13 )  max_pr_user_on_file
+    IF ( max_pr_user_on_file > max_pr_user )  THEN
+       IF ( myid == 0 )  THEN
+          PRINT*, '+++ read_parts_of_var_list: too many user profiles on res', &
+                       'tart data file'
+          PRINT*, '                   max_pr_user on file    = "', &
+                  max_pr_user_on_file, '"'
+          PRINT*, '                   max_pr_user from run   = "', &
+                  max_pr_user, '"'
+       ENDIF
+       CALL local_stop
+    ENDIF
+
+    READ ( 13 )  variable_chr
+    IF ( TRIM( variable_chr ) /= 'statistic_regions' )  THEN
+       PRINT*, '+++ read_var_list: statistic_regions not found in restart da', &
+                    'ta file'
+       CALL local_stop
+    ENDIF
+    READ ( 13 )  statistic_regions_on_file
+    IF ( statistic_regions_on_file > statistic_regions )  THEN
+       IF ( myid == 0 )  THEN
+          PRINT*, '+++ read_parts_of_var_list: too many statistic regions on', &
+                       ' restart data file'
+          PRINT*, '                   statistic regions on file    = "', &
+                  max_pr_user_on_file, '"'
+          PRINT*, '                   statistic regions from run   = "', &
+                  max_pr_user, '"'
+       ENDIF
+       CALL local_stop
+    ENDIF
+
+
+!
+!-- Now read and check some control parameters and skip the rest
+!-- The total domain of the pre-run must not be smaller than the subdomain
+!-- of the current run, because the mapping of data from the pre-run does
+!-- not work for this case.
+    WRITE (9,*) 'wpovl: begin reading variables'
+    CALL local_flush( 9 )
+    READ ( 13 )  variable_chr
+
+    DO  WHILE ( TRIM( variable_chr ) /= '*** end ***' )
+
+       SELECT CASE ( TRIM( variable_chr ) )
+
+          CASE ( 'hom' )
+             ALLOCATE( hom_on_file(0:nz+1,2,pr_palm+max_pr_user_on_file, &
+                       0:statistic_regions_on_file) )
+             READ ( 13 )  hom_on_file
+             hom = hom_on_file(:,:,1:pr_palm+max_pr_user,0:statistic_regions)
+             DEALLOCATE( hom_on_file )
+
+          CASE ( 'nx' )
+             READ ( 13 )  nx_on_file
+             IF ( nx_on_file < ( nxr - nxl ) )  THEN
+                PRINT*, '+++ read_parts_of_var_list: total domain along x on', &
+                             ' restart file is smaller than current subdomain'
+                PRINT*, '                            nx on file = ', nx_on_file
+                PRINT*, '                            nxr - nxl  = ', nxr - nxl
+                CALL local_stop
+             ENDIF
+
+          CASE ( 'ny' )
+             READ ( 13 )  ny_on_file
+             IF ( ny_on_file < ( nyn - nys ) )  THEN
+                PRINT*, '+++ read_parts_of_var_list: total domain along y on', &
+                             ' restart file is smaller than current subdomain'
+                PRINT*, '                            ny on file = ', ny_on_file
+                PRINT*, '                            nyn - nys  = ', nyn - nys
+                CALL local_stop
+             ENDIF
+
+
+          CASE DEFAULT
+
+             READ ( 13 )  idum
+
+       END SELECT
+
+       READ ( 13 )  variable_chr
+
+    ENDDO
+
+
+ END SUBROUTINE read_parts_of_var_list
+
+
+
  SUBROUTINE skip_var_list
 
@@ -550,4 +723,6 @@
 !------------------------------------------------------------------------------!
 
+    IMPLICIT NONE
+
     CHARACTER (LEN=10) ::  version_on_file
     CHARACTER (LEN=30) ::  variable_chr
@@ -556,14 +731,24 @@
 
 
+    WRITE (9,*) 'skipvl #1'
+    CALL local_flush( 9 )
     READ ( 13 )  version_on_file
 
-    READ ( 13 )  variable_chr
+    WRITE (9,*) 'skipvl before variable_chr'
+    CALL local_flush( 9 )
+    READ ( 13 )  variable_chr
+    WRITE (9,*) 'skipvl after variable_chr'
+    CALL local_flush( 9 )
 
     DO  WHILE ( TRIM( variable_chr ) /= '*** end ***' )
 
+    WRITE (9,*) 'skipvl chr = ', variable_chr
+    CALL local_flush( 9 )
        READ ( 13 )  idum
        READ ( 13 )  variable_chr
 
     ENDDO
+    WRITE (9,*) 'skipvl last'
+    CALL local_flush( 9 )
 
 
Index: /palm/trunk/SOURCE/user_interface.f90
===================================================================
--- /palm/trunk/SOURCE/user_interface.f90	(revision 146)
+++ /palm/trunk/SOURCE/user_interface.f90	(revision 147)
@@ -177,5 +177,7 @@
 
 
- SUBROUTINE user_read_restart_data
+ SUBROUTINE user_read_restart_data( nxlc, nxlf, nxl_on_file, nxrc, nxrf,  &
+                                    nxr_on_file, nync, nynf, nyn_on_file, &
+                                    nysc, nysf, nys_on_file, tmp_2d, tmp_3d )
 
 !------------------------------------------------------------------------------!
@@ -184,4 +186,8 @@
 ! ------------
 ! Reading restart data from file(s)
+! Subdomain index limits on file are given by nxl_on_file, etc.
+! Indices nxlc, etc. indicate the range of gridpoints to be mapped from the
+! subdomain on file (f) to the subdomain of the current PE (c). They have been
+! calculated in routine read_3d_binary.
 !------------------------------------------------------------------------------!
 
@@ -194,4 +200,15 @@
 
     CHARACTER (LEN=20) :: field_char
+
+    INTEGER ::  nxlc, nxlf, nxl_on_file, nxrc, nxrf, nxr_on_file, nync, nynf, &
+                nyn_on_file, nysc, nysf, nys_on_file
+
+    REAL, DIMENSION(nys_on_file-1:nyn_on_file+1,nxl_on_file-1:nxr_on_file+1) ::&
+          tmp_2d
+
+    REAL, DIMENSION(nzb:nzt+1,nys_on_file-1:nyn_on_file+1, &
+                    nxl_on_file-1:nxr_on_file+1) ::        &
+          tmp_3d
+
 !
 !-- Here the reading of user-defined restart data follows:
@@ -207,6 +224,11 @@
 !
 !             CASE ( 'u2_av' )
-!                ALLOCATE( u2_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+!                IF ( .NOT. ALLOCATED( u2_av )  THEN
+!                   ALLOCATE( u2_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+!                ENDIF
+!                READ ( 13 )  tmp_3d
 !                READ ( 13 )  u2_av
+!                u2_av(:,nysc-1:nync+1,nxlc-1:nxrc+1) = &
+!                                          tmp_3d(:,nysf-1:nynf+1,nxlf-1:nxrf+1)
 !
 !             CASE DEFAULT
