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 )
 
 !
