Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4897)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4898)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! MPI-IO for sky view factors implemented
+!
+! 4896 2021-03-03 16:10:18Z raasch
 ! typo in file appendix removed
 !
@@ -477,4 +480,5 @@
                pt_surface,                                                                         &
                read_svf,                                                                           &
+               restart_data_format_input,                                                          &
                restart_data_format_output,                                                         &
                rho_surface,                                                                        &
@@ -506,4 +510,5 @@
                nxr,                                                                                &
                nxrg,                                                                               &
+               nx_on_file,                                                                         &
                ny,                                                                                 &
                nyn,                                                                                &
@@ -511,4 +516,5 @@
                nys,                                                                                &
                nysg,                                                                               &
+               ny_on_file,                                                                         &
                nzb,                                                                                &
                nzt,                                                                                &
@@ -595,5 +601,11 @@
         ONLY:  rd_mpi_io_check_array,                                                              &
                rrd_mpi_io,                                                                         &
-               wrd_mpi_io
+               wrd_mpi_io,                                                                         &
+               rd_mpi_io_open,                                                                     &
+               rd_mpi_io_surface_filetypes,                                                        &
+               rrd_mpi_io_surface,                                                                 &
+               wrd_mpi_io_surface,                                                                 &
+               rd_mpi_io_close,                                                                    &
+               tgh
 
     USE statistics,                                                                                &
@@ -3424,35 +3436,36 @@
 
 !
-!-- Find all discretized apparent solar positions for radiation interaction.
-    IF ( radiation_interactions )  CALL radiation_presimulate_solar_pos
-
-!
-!-- If required, read or calculate and write out the SVF
-    IF ( radiation_interactions .AND. read_svf )  THEN
-!
-!--    Read sky-view factors and further required data from file
-       CALL radiation_read_svf()
-
-    ELSEIF ( radiation_interactions .AND. .NOT. read_svf )  THEN
-!
-!--    Calculate SFV and CSF
-       CALL radiation_calc_svf()
+!-- If required, read or calculate and write out the SVF.
+    IF ( radiation_interactions )  THEN
+
+!
+!--    Find all discretized apparent solar positions for radiation interaction.
+       CALL radiation_presimulate_solar_pos
+
+       IF ( read_svf )  THEN
+!
+!--       Read sky-view factors and further required data from file
+          CALL radiation_read_svf()
+
+       ELSE
+!
+!--       Calculate svf and csf.
+          CALL radiation_calc_svf()
+       ENDIF
+
+       IF ( write_svf )  THEN
+!
+!--       Write svf, csf svfsurf and csfsurf data to file.
+          CALL radiation_write_svf()
+       ENDIF
+
+!
+!--    Adjust radiative fluxes.
+!--    In case of urban and land surfaces, also call an initial interaction.
+       CALL radiation_interaction
+
     ENDIF
 
-    IF ( radiation_interactions .AND. write_svf )  THEN
-!
-!--    Write svf, csf svfsurf and csfsurf data to file
-       CALL radiation_write_svf()
-    ENDIF
-
-!
-!-- Adjust radiative fluxes. In case of urban and land surfaces, also call an initial interaction.
-    IF ( radiation_interactions )  THEN
-       CALL radiation_interaction
-    ENDIF
-
     IF ( debug_output )  CALL debug_message( 'radiation_init', 'end' )
-
-    RETURN ! @Todo: remove, I don't see what we need this for here
 
  END SUBROUTINE radiation_init
@@ -9450,142 +9463,436 @@
     CHARACTER(rad_version_len) ::  rad_version_field  !<
 
+    CHARACTER(LEN=64) ::  rad_version_in  !< rad version on file
+
     INTEGER(iwp) ::  i                      !<
+    INTEGER(iwp) ::  ncsfl_tot              !<
     INTEGER(iwp) ::  ndsidir_from_file = 0  !<
+    INTEGER(iwp) ::  ndsidir_tot            !<
     INTEGER(iwp) ::  npcbl_from_file   = 0  !<
+    INTEGER(iwp) ::  npcbl_tot              !<
     INTEGER(iwp) ::  nsurfl_from_file  = 0  !<
+    INTEGER(iwp) ::  nsurfl_tot             !<
+    INTEGER(iwp) ::  nsvfl_tot              !< total (sum across all PEs) counter for the different IO variables
     INTEGER(iwp) ::  nmrtbl_from_file  = 0  !<
-
-
-    CALL location_message( 'reading view factors for radiation interaction', 'start' )
-
-    DO  i = 0, io_blocks-1
-       IF ( i == io_group )  THEN
-
-!
-!--       numprocs_previous_run is only known in case of reading restart data. If a new initial run
-!--       which reads svf data is started the following query will be skipped.
-          IF ( initializing_actions == 'read_restart_data' )  THEN
-
-             IF ( numprocs_previous_run /= numprocs )  THEN
-                WRITE( message_string, * ) 'A different number of processors between the run ',    &
-                                           'that has written the svf data and the one that ',      &
-                                           'will read it is not allowed'
-                CALL message( 'check_open', 'PA0491', 1, 2, 0, 6, 0 )
+    INTEGER(iwp) ::  nmrtbl_tot             !<
+    INTEGER(iwp) ::  nmrtf_tot              !<
+
+    INTEGER(iwp), DIMENSION(4) ::  global_sum  !<
+    INTEGER(iwp), DIMENSION(4) ::  local_sum   !< variables to compute total counter
+
+    INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) ::  global_end    !< global end index
+    INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) ::  global_start  !< global start index
+    INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) ::  end_index     !< local end index
+    INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) ::  start_index   !< local start index
+
+    LOGICAL ::  data_to_read  !< flag indicating if data is available for current variable
+
+    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  tmp  !<
+
+
+    IF ( TRIM( restart_data_format_input ) == 'fortran_binary' )  THEN
+
+       CALL location_message( 'reading sky view factors for radiation interaction ' //             &
+                              '(Fortran binary format)', 'start' )
+
+       DO  i = 0, io_blocks-1
+          IF ( i == io_group )  THEN
+
+!
+!--          numprocs_previous_run is only known in case of reading restart data. If a new initial
+!--          run which reads svf data is started the following query will be skipped.
+             IF ( initializing_actions == 'read_restart_data' )  THEN
+
+                IF ( numprocs_previous_run /= numprocs )  THEN
+                   WRITE( message_string, * ) 'A different number of processors between the run ', &
+                                              'that has written the svf data and the one that ',   &
+                                              'will read it is not allowed'
+                   CALL message( 'check_open', 'PA0491', 1, 2, 0, 6, 0 )
+                ENDIF
+
              ENDIF
 
+!
+!--          Open binary file
+             CALL check_open( 88 )
+
+!
+!--          Read and check version
+             READ ( 88 ) rad_version_field
+             IF ( TRIM( rad_version_field ) /= TRIM( rad_version ) )  THEN
+                 WRITE( message_string, * ) 'Version of binary SVF file "',                        &
+                                            TRIM( rad_version_field ), '" does not match ',        &
+                                            'the version of model "', TRIM( rad_version ), '"'
+                 CALL message( 'radiation_read_svf', 'PA0482', 1, 2, 0, 6, 0 )
+             ENDIF
+
+!
+!--          Read nsvfl, ncsfl, nsurfl, nmrtf
+             READ ( 88 ) nsvfl, ncsfl, nsurfl_from_file, npcbl_from_file, ndsidir_from_file,       &
+                         nmrtbl_from_file, nmrtf
+
+             IF ( nsvfl < 0  .OR.  ncsfl < 0 )  THEN
+                 WRITE( message_string, * ) 'Wrong number of SVF or CSF'
+                 CALL message( 'radiation_read_svf', 'PA0483', 1, 2, 0, 6, 0 )
+             ELSE
+                 WRITE(debug_string,*) 'Number of SVF, CSF, and nsurfl to read', nsvfl, ncsfl,     &
+                                       nsurfl_from_file
+                 IF ( debug_output )  CALL debug_message( debug_string, 'info' )
+             ENDIF
+
+             IF ( nsurfl_from_file /= nsurfl )  THEN
+                 WRITE( message_string, * ) 'nsurfl from SVF file does not match calculated ',     &
+                                            'nsurfl from radiation_interaction_init'
+                 CALL message( 'radiation_read_svf', 'PA0490', 1, 2, 0, 6, 0 )
+             ENDIF
+
+             IF ( npcbl_from_file /= npcbl )  THEN
+                 WRITE( message_string, * ) 'npcbl from SVF file does not match calculated npcbl', &
+                                            ' from radiation_interaction_init'
+                 CALL message( 'radiation_read_svf', 'PA0493', 1, 2, 0, 6, 0 )
+             ENDIF
+
+             IF ( ndsidir_from_file /= ndsidir )  THEN
+                 WRITE( message_string, * ) 'ndsidir from SVF file does not match calculated ',    &
+                                            'ndsidir from radiation_presimulate_solar_pos'
+                 CALL message( 'radiation_read_svf', 'PA0494', 1, 2, 0, 6, 0 )
+             ENDIF
+
+             IF ( nmrtbl_from_file /= nmrtbl )  THEN
+                 WRITE( message_string, * ) 'nmrtbl from SVF file does not match calculated ',     &
+                                            'nmrtbl from radiation_interaction_init'
+                 CALL message( 'radiation_read_svf', 'PA0494', 1, 2, 0, 6, 0 )
+             ELSE
+                 WRITE( debug_string, * ) 'Number of nmrtf to read ', nmrtf
+                 IF ( debug_output )  CALL debug_message( debug_string, 'info' )
+             ENDIF
+
+!
+!--          Arrays skyvf, skyvft, dsitrans and dsitransc are allready allocated in
+!--          radiation_interaction_init and radiation_presimulate_solar_pos
+             IF ( nsurfl > 0 )  THEN
+                READ( 88 ) skyvf
+                READ( 88 ) skyvft
+                READ( 88 ) dsitrans
+             ENDIF
+
+             IF ( plant_canopy  .AND.  npcbl > 0 )  THEN
+                READ( 88 )  dsitransc
+             ENDIF
+
+!
+!--          The allocation of svf, svfsurf, csf, csfsurf, mrtf, mrtft, and mrtfsurf happens in
+!--          routine radiation_calc_svf which is not called if the program enters
+!--          radiation_read_svf. Therefore these arrays have to be allocated in the following.
+             IF ( nsvfl > 0 )  THEN
+                ALLOCATE( svf(ndsvf,nsvfl) )
+                ALLOCATE( svfsurf(idsvf,nsvfl) )
+                READ( 88 ) svf
+                READ( 88 ) svfsurf
+             ENDIF
+
+             IF ( plant_canopy  .AND.  ncsfl > 0 )  THEN
+                ALLOCATE( csf(ndcsf,ncsfl) )
+                ALLOCATE( csfsurf(idcsf,ncsfl) )
+                READ( 88 ) csf
+                READ( 88 ) csfsurf
+             ENDIF
+
+             IF ( nmrtbl > 0 )  THEN
+                READ( 88 ) mrtsky
+                READ( 88 ) mrtskyt
+                READ( 88 ) mrtdsit
+             ENDIF
+
+             IF ( nmrtf > 0 )  THEN
+                ALLOCATE( mrtf(nmrtf) )
+                ALLOCATE( mrtft(nmrtf) )
+                ALLOCATE( mrtfsurf(2,nmrtf) )
+                READ( 88 ) mrtf
+                READ( 88 ) mrtft
+                READ( 88 ) mrtfsurf
+             ENDIF
+
+!
+!--          Close binary file
+             CALL close_file( 88 )
+
           ENDIF
-
-!
-!--       Open binary file
-          CALL check_open( 88 )
-
-!
-!--       Read and check version
-          READ ( 88 ) rad_version_field
-          IF ( TRIM( rad_version_field ) /= TRIM( rad_version ) )  THEN
-              WRITE( message_string, * ) 'Version of binary SVF file "',                           &
-                                         TRIM( rad_version_field ), '" does not match ',           &
-                                         'the version of model "', TRIM( rad_version ), '"'
-              CALL message( 'radiation_read_svf', 'PA0482', 1, 2, 0, 6, 0 )
+#if defined( __parallel )
+          CALL MPI_BARRIER( comm2d, ierr )
+#endif
+       ENDDO
+
+       CALL location_message( 'reading sky view factors for radiation interaction ' //             &
+                              '(Fortran binary format)', 'finished' )
+
+
+    ELSEIF ( restart_data_format_input(1:3) == 'mpi' )  THEN
+
+!
+!--    In case of MPI-IO data is treated like surface data and the respective routines from
+!--    restart_data_mpi_io_mod are used for reading. Contrary to restart data, the virtual processor
+!--    grid can not be changed between writing and reading svf data.
+       CALL location_message( 'reading sky view factors for radiation interaction ' //             &
+                              '(MPI-IO format)', 'start' )
+!
+!--    Open MPI-IO svf file for reading the global data.
+       CALL rd_mpi_io_open( 'READ', 'SVFIN' // TRIM( coupling_char ),                              &
+                            open_for_global_io_only = .TRUE. )
+!
+!--    Check general header.
+       IF ( tgh%pes_along_x /= npex  .OR.  tgh%pes_along_y /= npey )  THEN
+!
+!--       Force re-calculation of svfs.
+          read_svf = .FALSE.
+          WRITE( message_string, '(A,I7,A,I7,A,I7,A,I7,A)' )                                       &
+              'virtual PE grid has changed between previous and current run &npex_prev = ',        &
+              tgh%pes_along_x, ' npey_prev = ', tgh%pes_along_y, ' npex_new = ', npex,             &
+              ' npey_new = ', npey, '&svf wil be re-calculated'
+          CALL message( 'radiation_read_svf', 'PA0517', 0, 0, 0, 6, 0 )
+          RETURN
+       ENDIF
+!
+!--    Read global variables
+       CALL rrd_mpi_io( 'rad_version', rad_version_in )
+       CALL rrd_mpi_io( 'nsvfl', nsvfl_tot )
+       CALL rrd_mpi_io( 'ncsfl', ncsfl_tot )
+       CALL rrd_mpi_io( 'nsurfl', nsurfl_from_file )
+       CALL rrd_mpi_io( 'npcbl', npcbl_from_file )
+       CALL rrd_mpi_io( 'ndsidir', ndsidir_from_file )
+       CALL rrd_mpi_io( 'nmrtbl', nmrtbl_from_file )
+       CALL rrd_mpi_io( 'nmrtf', nmrtf_tot )
+
+       CALL rd_mpi_io_close
+!
+!--    Compute global values of local counters.
+       local_sum(1) = nsurfl
+       local_sum(2) = npcbl
+       local_sum(3) = ndsidir
+       local_sum(4) = nmrtbl
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( local_sum, global_sum, SIZE(local_sum), MPI_INTEGER, MPI_SUM, comm2d,   &
+                           ierr)
+#else
+       global_sum = local_sum
+#endif
+       nsurfl_tot  = global_sum(1)
+       npcbl_tot   = global_sum(2)
+       ndsidir_tot = global_sum(3)
+       nmrtbl_tot  = global_sum(4)
+!
+!--    Check for errors.
+       nx_on_file = tgh%total_nx-1
+       ny_on_file = tgh%total_ny-1
+
+       IF ( nx_on_file /= nx  .OR.  ny_on_file /= ny )  THEN
+          WRITE( message_string, '(A,4(A,I7))' )                                                   &
+               'total number of grid points along x and y in file SVFIN do not match current run', &
+               '&nx_on_file = ', nx_on_file, ' ny_on_file = ', ny_on_file, ' nx = ', nx, ' ny = ', &
+               ny
+          CALL message( 'radiation_read_svf', 'PA0518', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( TRIM( rad_version_in ) /= TRIM( rad_version ) )  THEN
+          WRITE( message_string, * ) 'Version of binary SVF file "', TRIM( rad_version_field ),    &
+                                     '" does not match the version of model "',                    &
+                                     TRIM( rad_version ), '"'
+          CALL message( 'radiation_read_svf', 'PA0482', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( nsvfl_tot < 0  .OR.  ncsfl_tot < 0 )  THEN
+          WRITE( message_string, * ) 'Wrong number of SVF or CSF'
+          CALL message( 'radiation_read_svf', 'PA0483', 1, 2, 0, 6, 0 )
+       ELSE
+          WRITE(debug_string,*) 'Number of SVF, CSF, and nsurfl to read', nsvfl_tot, ncsfl_tot,  &
+                                          nsurfl_tot
+          IF ( debug_output )  CALL debug_message( debug_string, 'info' )
+       ENDIF
+
+       IF ( nsurfl_from_file /= nsurfl_tot )  THEN
+          WRITE( message_string, * ) 'nsurfl from SVF file does not match calculated ',        &
+             'nsurfl from radiation_interaction_init'
+          CALL message( 'radiation_read_svf', 'PA0490', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( npcbl_from_file /= npcbl_tot )  THEN
+          WRITE( message_string, * ) 'npcbl from SVF file does not match calculated npcbl ',   &
+             'from radiation_interaction_init'
+          CALL message( 'radiation_read_svf', 'PA0493', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( ndsidir_from_file /= ndsidir_tot )  THEN
+          WRITE( message_string, * ) 'ndsidir from SVF file does not match calculated ',       &
+             'ndsidir from radiation_presimulate_solar_pos'
+          CALL message( 'radiation_read_svf', 'PA0494', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( nmrtbl_from_file /= nmrtbl_tot )  THEN
+          WRITE( message_string, * ) 'nmrtbl from SVF file does not match calculated nmrtbl ', &
+             'from radiation_interaction_init'
+          CALL message( 'radiation_read_svf', 'PA0494', 1, 2, 0, 6, 0 )
+       ELSE
+          WRITE( debug_string, * ) 'Number of nmrtf to read ', nmrtf_tot
+          IF ( debug_output )  CALL debug_message( debug_string, 'info' )
+       ENDIF
+
+!
+!--    Open MPI-IO SVF file for read local data.
+       CALL rd_mpi_io_open( 'READ', 'SVFIN' // TRIM( coupling_char ) )
+
+       IF ( nsurfl_tot > 0 )  THEN
+!
+!--       Read global indices.
+          CALL rrd_mpi_io( 'nsurfl_global_start', global_start )
+          CALL rrd_mpi_io( 'nsurfl_global_end', global_end )
+!
+!--       Set file types of variables and compute local indices.
+          CALL rd_mpi_io_surface_filetypes( start_index, end_index, data_to_read, global_start,    &
+                                            global_end )
+          nsurfl = end_index(nyn,nxr)
+
+          IF ( data_to_read )  THEN
+
+             CALL rrd_mpi_io_surface( 'skyvf', skyvf )
+             CALL rrd_mpi_io_surface( 'skyvft', skyvft )
+!
+!--          To avoid another overlay of rrd_mpi_io_surface, dsitrans is read as REAL tmp array.
+!--          The order of dimensions of dsitrans is different to the order expected by
+!--          rrd_mpi_io_surface. Therefor a tranpose of tmp is required.
+             ALLOCATE( tmp(SIZE(dsitrans,2),SIZE(dsitrans,1)) )
+             CALL rrd_mpi_io_surface( 'dsitrans', tmp )
+             dsitrans = TRANSPOSE( tmp )
+             DEALLOCATE( tmp )
+
           ENDIF
 
-!
-!--       Read nsvfl, ncsfl, nsurfl, nmrtf
-          READ ( 88 ) nsvfl, ncsfl, nsurfl_from_file, npcbl_from_file, ndsidir_from_file,          &
-                      nmrtbl_from_file, nmrtf
-
-          IF ( nsvfl < 0  .OR.  ncsfl < 0 )  THEN
-              WRITE( message_string, * ) 'Wrong number of SVF or CSF'
-              CALL message( 'radiation_read_svf', 'PA0483', 1, 2, 0, 6, 0 )
-          ELSE
-              WRITE(debug_string,*) 'Number of SVF, CSF, and nsurfl to read', nsvfl, ncsfl,        &
-                                    nsurfl_from_file
-              IF ( debug_output )  CALL debug_message( debug_string, 'info' )
+       ENDIF
+
+       IF ( npcbl_tot > 0 )  THEN
+
+          CALL rrd_mpi_io( 'npcbl_global_start', global_start )
+          CALL rrd_mpi_io( 'npcbl_global_end', global_end )
+
+          CALL rd_mpi_io_surface_filetypes( start_index, end_index, data_to_read, global_start,    &
+                                            global_end )
+          npcbl = end_index(nyn,nxr)
+
+          IF ( data_to_read )  THEN
+
+             ALLOCATE( tmp(SIZE(dsitransc,2),SIZE(dsitransc,1)) )
+             CALL rrd_mpi_io_surface( 'dsitransc', tmp )
+             dsitransc = TRANSPOSE( tmp )
+             DEALLOCATE( tmp )
+
           ENDIF
 
-          IF ( nsurfl_from_file /= nsurfl )  THEN
-              WRITE( message_string, * ) 'nsurfl from SVF file does not match calculated ',        &
-                                         'nsurfl from radiation_interaction_init'
-              CALL message( 'radiation_read_svf', 'PA0490', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( nsvfl_tot > 0 )  THEN
+
+          CALL rrd_mpi_io( 'nsvfl_global_start', global_start )
+          CALL rrd_mpi_io( 'nsvfl_global_end', global_end )
+
+          CALL rd_mpi_io_surface_filetypes( start_index, end_index, data_to_read, global_start,    &
+                                            global_end )
+          nsvfl = end_index(nyn,nxr)
+
+          IF ( .NOT. ALLOCATED( svf )     )  ALLOCATE( svf(ndsvf,nsvfl) )
+          IF ( .NOT. ALLOCATED( svfsurf ) )  ALLOCATE( svfsurf(idsvf,nsvfl) )
+
+          IF ( data_to_read )  THEN
+
+             CALL rrd_mpi_io_surface( 'svf', svf )
+             ALLOCATE( tmp(SIZE(svfsurf,1),SIZE(svfsurf,2)) )
+             CALL rrd_mpi_io_surface( 'svfsurf', tmp )
+             svfsurf = tmp
+             DEALLOCATE( tmp )
+
           ENDIF
 
-          IF ( npcbl_from_file /= npcbl )  THEN
-              WRITE( message_string, * ) 'npcbl from SVF file does not match calculated npcbl ',   &
-                                         'from radiation_interaction_init'
-              CALL message( 'radiation_read_svf', 'PA0493', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( plant_canopy  )  THEN
+
+          IF ( ncsfl_tot > 0 )  THEN
+
+             CALL rrd_mpi_io( 'ncsfl_global_start', global_start )
+             CALL rrd_mpi_io( 'ncsfl_global_end', global_end )
+
+             CALL rd_mpi_io_surface_filetypes( start_index, end_index, data_to_read, global_start, &
+                                               global_end )
+             ncsfl = end_index(nyn,nxr)
+
+             IF ( .NOT. ALLOCATED( csf )     )  ALLOCATE( csf(ndcsf,ncsfl) )
+             IF ( .NOT. ALLOCATED( csfsurf ) )  ALLOCATE( csfsurf(idcsf,ncsfl) )
+
+             IF ( data_to_read )  THEN
+
+                CALL rrd_mpi_io_surface( 'csf', csf )
+                ALLOCATE( tmp(SIZE(csfsurf,1),SIZE(csfsurf,2)) )
+                CALL rrd_mpi_io_surface( 'csfsurf', tmp )
+                csfsurf = tmp
+                DEALLOCATE( tmp )
+
+             ENDIF
+
           ENDIF
 
-          IF ( ndsidir_from_file /= ndsidir )  THEN
-              WRITE( message_string, * ) 'ndsidir from SVF file does not match calculated ',       &
-                                         'ndsidir from radiation_presimulate_solar_pos'
-              CALL message( 'radiation_read_svf', 'PA0494', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( nmrtbl_tot > 0 )  THEN
+
+          CALL rrd_mpi_io( 'nmrtbl_global_start', global_start )
+          CALL rrd_mpi_io( 'nmrtbl_global_end', global_end )
+
+          CALL rd_mpi_io_surface_filetypes( start_index, end_index, data_to_read, global_start,    &
+                                            global_end )
+          nmrtbl = end_index(nyn,nxr)
+
+          IF ( data_to_read )  THEN
+
+             CALL rrd_mpi_io_surface( 'mrtsky', mrtsky )
+             CALL rrd_mpi_io_surface( 'mrtskyt', mrtskyt )
+             ALLOCATE( tmp(SIZE(mrtdsit,2),SIZE(mrtdsit,1)) )
+             CALL rrd_mpi_io_surface( 'mrtdsit', tmp )
+             mrtdsit = TRANSPOSE( tmp )
+             DEALLOCATE( tmp )
+
           ENDIF
-          IF ( nmrtbl_from_file /= nmrtbl )  THEN
-              WRITE( message_string, * ) 'nmrtbl from SVF file does not match calculated nmrtbl ', &
-                                         'from radiation_interaction_init'
-              CALL message( 'radiation_read_svf', 'PA0494', 1, 2, 0, 6, 0 )
-          ELSE
-              WRITE( debug_string, * ) 'Number of nmrtf to read ', nmrtf
-              IF ( debug_output )  CALL debug_message( debug_string, 'info' )
+
+       ENDIF
+
+       IF ( nmrtf_tot > 0 )  THEN
+
+          CALL rrd_mpi_io( 'nmrtf_global_start', global_start )
+          CALL rrd_mpi_io( 'nmrtf_global_end', global_end )
+
+          CALL rd_mpi_io_surface_filetypes( start_index, end_index, data_to_read, global_start,    &
+                                            global_end )
+          nmrtf = end_index(nyn,nxr)
+
+          IF ( .NOT. ALLOCATED( mrtf )     )  ALLOCATE( mrtf(nmrtf) )
+          IF ( .NOT. ALLOCATED( mrtft )    )  ALLOCATE( mrtft(nmrtf) )
+          IF ( .NOT. ALLOCATED( mrtfsurf ) )  ALLOCATE( mrtfsurf(2,nmrtf) )
+
+          IF ( data_to_read )  THEN
+
+             CALL rrd_mpi_io_surface( 'mrtf', mrtf )
+             CALL rrd_mpi_io_surface( 'mrtft', mrtft )
+             ALLOCATE( tmp(SIZE(mrtfsurf,1),SIZE(mrtfsurf,2)) )
+             CALL rrd_mpi_io_surface( 'mrtfsurf', tmp )
+             mrtfsurf = tmp
+             DEALLOCATE( tmp )
+
           ENDIF
 
-!
-!--       Arrays skyvf, skyvft, dsitrans and dsitransc are allready allocated in
-!--       radiation_interaction_init and radiation_presimulate_solar_pos
-          IF ( nsurfl > 0 )  THEN
-             READ( 88 ) skyvf
-             READ( 88 ) skyvft
-             READ( 88 ) dsitrans
-          ENDIF
-
-          IF ( plant_canopy  .AND.  npcbl > 0 )  THEN
-             READ( 88 )  dsitransc
-          ENDIF
-
-!
-!--       The allocation of svf, svfsurf, csf, csfsurf, mrtf, mrtft, and mrtfsurf happens in
-!--       routine radiation_calc_svf which is not called if the program enters radiation_read_svf.
-!--       Therefore these arrays have to be allocated in the following.
-          IF ( nsvfl > 0 )  THEN
-             ALLOCATE( svf(ndsvf,nsvfl) )
-             ALLOCATE( svfsurf(idsvf,nsvfl) )
-             READ( 88 ) svf
-             READ( 88 ) svfsurf
-          ENDIF
-
-          IF ( plant_canopy  .AND.  ncsfl > 0 )  THEN
-             ALLOCATE( csf(ndcsf,ncsfl) )
-             ALLOCATE( csfsurf(idcsf,ncsfl) )
-             READ( 88 ) csf
-             READ( 88 ) csfsurf
-          ENDIF
-
-          IF ( nmrtbl > 0 )  THEN
-             READ( 88 ) mrtsky
-             READ( 88 ) mrtskyt
-             READ( 88 ) mrtdsit
-          ENDIF
-
-          IF ( nmrtf > 0 )  THEN
-             ALLOCATE( mrtf(nmrtf) )
-             ALLOCATE( mrtft(nmrtf) )
-             ALLOCATE( mrtfsurf(2,nmrtf) )
-             READ( 88 ) mrtf
-             READ( 88 ) mrtft
-             READ( 88 ) mrtfsurf
-          ENDIF
-
-!
-!--       Close binary file
-          CALL close_file( 88 )
-
-       ENDIF
-#if defined( __parallel )
-       CALL MPI_BARRIER( comm2d, ierr )
-#endif
-    ENDDO
-
-    CALL location_message( 'reading view factors for radiation interaction', 'finished' )
-
+       ENDIF
+
+       CALL rd_mpi_io_close
+
+       CALL location_message( 'reading sky view factors for radiation interaction ' //             &
+                              '(MPI-IO format)', 'finished' )
+
+    ENDIF
 
  END SUBROUTINE radiation_read_svf
@@ -9604,55 +9911,444 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i  !<
-
-
-    CALL location_message( 'writing view factors for radiation interaction', 'start' )
-
-    DO  i = 0, io_blocks-1
-       IF ( i == io_group )  THEN
-!
-!--       Open binary file
-          CALL check_open( 89 )
-
-          WRITE( 89 ) rad_version
-          WRITE( 89 ) nsvfl, ncsfl, nsurfl, npcbl, ndsidir, nmrtbl, nmrtf
-          IF ( nsurfl > 0 )  THEN
-             WRITE( 89 ) skyvf
-             WRITE( 89 ) skyvft
-             WRITE( 89 ) dsitrans
+    INTEGER(iwp), PARAMETER ::  max_i4_value  = 2147483647  !< maximum positive INTEGER(4) value, 2**31 - 1
+
+    INTEGER(iwp) ::  i            !<
+    INTEGER(iwp) ::  ic           !<
+    INTEGER(iwp) ::  ierr         !<
+    INTEGER(iwp) ::  ind          !<
+    INTEGER(iwp) ::  ipcgb        !<
+    INTEGER(iwp) ::  isurf        !<
+    INTEGER(iwp) ::  j            !<
+    INTEGER(iwp) ::  jc           !<
+    INTEGER(iwp) ::  ncsfl_tot    !<
+    INTEGER(iwp) ::  ndsidir_tot  !<
+    INTEGER(iwp) ::  nmrtbl_tot   !<
+    INTEGER(iwp) ::  nmrtf_tot    !<
+    INTEGER(iwp) ::  npcbl_tot    !<
+    INTEGER(iwp) ::  nsurfl_tot   !<
+    INTEGER(iwp) ::  nsvfl_tot    !< total (sum over all PEs) counter for the different IO variables
+
+    INTEGER(idp), DIMENSION(7) ::  local_sum   !< variables to compute total counter
+    INTEGER(idp), DIMENSION(7) ::  global_sum  !< INTEGER(idp) to allow check, if total number of values > 2G
+
+    INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) ::  end_index     !< local end index
+    INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) ::  global_end    !< global end index
+    INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) ::  global_start  !< global start index
+    INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) ::  lo_no         !< local number of values
+    INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) ::  start_index   !< local start index
+
+    LOGICAL ::  data_to_write  !<  flag indicating if data is available for writing
+
+    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  tmp  !<
+
+
+    IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
+
+       CALL location_message( 'writing sky view factors for radiation interaction ' //             &
+                              '(Fortran binary format)', 'start' )
+
+       DO  i = 0, io_blocks-1
+          IF ( i == io_group )  THEN
+!
+!--          Open binary file
+             CALL check_open( 89 )
+
+             WRITE( 89 ) rad_version
+             WRITE( 89 ) nsvfl, ncsfl, nsurfl, npcbl, ndsidir, nmrtbl, nmrtf
+             IF ( nsurfl > 0 )  THEN
+                WRITE( 89 ) skyvf
+                WRITE( 89 ) skyvft
+                WRITE( 89 ) dsitrans
+             ENDIF
+             IF ( npcbl > 0 ) THEN
+                WRITE( 89 ) dsitransc
+             ENDIF
+             IF ( nsvfl > 0 ) THEN
+                WRITE( 89 ) svf
+                WRITE( 89 ) svfsurf
+             ENDIF
+             IF ( plant_canopy  .AND.  ncsfl > 0 )  THEN
+                 WRITE( 89 ) csf
+                 WRITE( 89 ) csfsurf
+             ENDIF
+             IF ( nmrtbl > 0 )  THEN
+                WRITE( 89 ) mrtsky
+                WRITE( 89 ) mrtskyt
+                WRITE( 89 ) mrtdsit
+             ENDIF
+             IF ( nmrtf > 0 )  THEN
+                 WRITE( 89 ) mrtf
+                 WRITE( 89 ) mrtft
+                 WRITE( 89 ) mrtfsurf
+             ENDIF
+!
+!--          Close binary file
+             CALL close_file( 89 )
+
           ENDIF
-          IF ( npcbl > 0 ) THEN
-             WRITE( 89 ) dsitransc
+#if defined( __parallel )
+          CALL MPI_BARRIER( comm2d, ierr )
+#endif
+       ENDDO
+
+       CALL location_message( 'writing sky view factors for radiation interaction ' //             &
+                              '(Fortran binary format)', 'fnished' )
+
+    ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
+!
+!--    Sky view factor data is treated like surface data and the respective routines from
+!--    restart_data_mpi_io_mod are used for writing. In contrary to restart data, the virtual
+!--    processor grid can not be changed between writing and reading svf data.
+       CALL location_message( 'writing sky view factors for radiation interaction ' //             &
+                              '(MPI-IO format)', 'start' )
+
+!
+!--    Open MPI-IO svf file
+       CALL rd_mpi_io_open( 'write', 'SVFOUT' // TRIM( coupling_char ) )
+!
+!--    Write global variables.
+       CALL wrd_mpi_io( 'rad_version', rad_version )
+!
+!--    Sum local number of skyview factor values on all PEs.
+       local_sum(1) = nsvfl
+       local_sum(2) = ncsfl
+       local_sum(3) = nsurfl
+       local_sum(4) = npcbl
+       local_sum(5) = ndsidir
+       local_sum(6) = nmrtbl
+       local_sum(7) = nmrtf
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( local_sum, global_sum, SIZE( local_sum ), MPI_INTEGER8, MPI_SUM,        &
+                           comm2d, ierr)
+#else
+       global_sum = local_sum
+#endif
+!
+!--    Check, if total number of respective skyview values do not exceed 2**31-1
+       DO  i = 1, 7
+          IF ( global_sum(i) > max_i4_value )  THEN
+             WRITE( message_string, '(A,I1,A)' ) 'number of sky view factor values (', i,          &
+                                                 ') > 2**31-1'
+             CALL message( 'radiation_write_svf', 'PA0733', 1, 2, 0, 6, 0 )
           ENDIF
-          IF ( nsvfl > 0 ) THEN
-             WRITE( 89 ) svf
-             WRITE( 89 ) svfsurf
+       ENDDO
+
+       nsvfl_tot   = global_sum(1)
+       ncsfl_tot   = global_sum(2)
+       nsurfl_tot  = global_sum(3)
+       npcbl_tot   = global_sum(4)
+       ndsidir_tot = global_sum(5)
+       nmrtbl_tot  = global_sum(6)
+       nmrtf_tot   = global_sum(7)
+!
+!--    Write total counters in header section of MPI-IO file.
+       CALL wrd_mpi_io( 'nsvfl', nsvfl_tot )
+       CALL wrd_mpi_io( 'ncsfl', ncsfl_tot )
+       CALL wrd_mpi_io( 'nsurfl', nsurfl_tot )
+       CALL wrd_mpi_io( 'npcbl', npcbl_tot )
+       CALL wrd_mpi_io( 'ndsidir', ndsidir_tot )
+       CALL wrd_mpi_io( 'nmrtbl', nmrtbl_tot )
+       CALL wrd_mpi_io( 'nmrtf', nmrtf_tot )
+!
+!--    Write local data.
+!--    All svf values are treated as surface values and use the respective routines from
+!--    restart_data_mpi_io_mod.
+       IF ( nsurfl > 0 )  THEN
+
+          lo_no = 0
+!
+!--       Count surface values on individual grid cells
+          DO  isurf = 1, nsurfl
+             jc = surfl(iy, isurf)
+             ic = surfl(ix, isurf)
+             lo_no(jc,ic) = lo_no(jc,ic) + 1
+          ENDDO
+!
+!--       Create local index array similar to surface routines.
+          ind = 1
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                start_index(j,i) = ind
+                end_index(j,i)   = start_index(j,i) + lo_no(j,i) - 1
+                ind = ind + lo_no(j,i)
+             ENDDO
+          ENDDO
+
+       ELSE
+
+          start_index = 1
+          end_index = 0
+
+       ENDIF
+!
+!--    Each PE has to call the next block, therefore nsurfl_tot is used.
+!--    This is required to use MPI_FILE_WRITE_All for writing.
+       IF ( nsurfl_tot > 0 )  THEN
+!
+!--       Set file types of variables for this block and compute global indices.
+          CALL rd_mpi_io_surface_filetypes( start_index, end_index, data_to_write, global_start,   &
+                                            global_end )
+
+          CALL wrd_mpi_io( 'nsurfl_global_start', global_start )
+          CALL wrd_mpi_io( 'nsurfl_global_end', global_end )
+
+          IF ( data_to_write )  THEN
+
+             CALL wrd_mpi_io_surface( 'skyvf', skyvf )
+             CALL wrd_mpi_io_surface( 'skyvft', skyvft )
+!
+!--          To avoid another overlay of rrd_mpi_io_surface, dsitrans is written as REAL tmp array.
+!--          The order of dimensions of dsitrans is different to the order expected by
+!--          rrd_mpi_io_surface. Therefor a tranpose of tmp is required.
+             ALLOCATE( tmp(SIZE(dsitrans,2),SIZE(dsitrans,1)) )
+             tmp = TRANSPOSE( dsitrans )
+             CALL wrd_mpi_io_surface( 'dsitrans', tmp )
+             DEALLOCATE( tmp )
+
           ENDIF
-          IF ( plant_canopy  .AND.  ncsfl > 0 )  THEN
-              WRITE( 89 ) csf
-              WRITE( 89 ) csfsurf
+
+       ENDIF
+
+       IF ( npcbl > 0 )  THEN
+
+          lo_no = 0
+          DO  isurf = 1, npcbl
+             jc = pcbl(iy,isurf)
+             ic = pcbl(ix,isurf)
+             lo_no(jc,ic) = lo_no(jc,ic) + 1
+          ENDDO
+
+          ind = 1
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                start_index(j,i) = ind
+                end_index(j,i)   = start_index(j,i) + lo_no(j,i) - 1
+                ind = ind+lo_no(j,i)
+             ENDDO
+          ENDDO
+
+       ELSE
+
+          start_index = 1
+          end_index = 0
+
+       ENDIF
+
+       IF ( npcbl_tot > 0 )  THEN
+
+          CALL rd_mpi_io_surface_filetypes( start_index, end_index, data_to_write, global_start,   &
+                                            global_end )
+          CALL wrd_mpi_io( 'npcbl_global_start', global_start )
+          CALL wrd_mpi_io( 'npcbl_global_end', global_end )
+
+          IF ( data_to_write )  THEN
+
+             ALLOCATE( tmp(SIZE(dsitransc,2),SIZE(dsitransc,1)) )
+             tmp = TRANSPOSE( dsitransc )
+             CALL wrd_mpi_io_surface( 'dsitransc', tmp )
+             DEALLOCATE( tmp )
           ENDIF
-          IF ( nmrtbl > 0 )  THEN
-             WRITE( 89 ) mrtsky
-             WRITE( 89 ) mrtskyt
-             WRITE( 89 ) mrtdsit
+       ENDIF
+
+       IF ( nsvfl > 0 )  THEN
+
+          lo_no = 0
+          DO  j = 1, SIZE( svfsurf, 2 )
+             isurf = svfsurf(1,j)
+             jc = surfl(iy,isurf)
+             ic = surfl(ix,isurf)
+             lo_no(jc,ic) = lo_no(jc,ic) + 1
+          ENDDO
+
+          ind = 1
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                start_index(j,i) = ind
+                end_index(j,i)   = start_index(j,i) + lo_no(j,i) - 1
+                ind = ind+lo_no(j,i)
+             ENDDO
+          ENDDO
+
+       ELSE
+
+          start_index = 1
+          end_index = 0
+       ENDIF
+
+       IF ( nsvfl_tot > 0 )  THEN
+
+          CALL rd_mpi_io_surface_filetypes( start_index, end_index, data_to_write, global_start,   &
+                                            global_end )
+          CALL wrd_mpi_io( 'nsvfl_global_start', global_start )
+          CALL wrd_mpi_io( 'nsvfl_global_end', global_end )
+
+          IF ( data_to_write )  THEN
+
+             CALL wrd_mpi_io_surface( 'svf', svf )
+             ALLOCATE( tmp(SIZE(svfsurf,1),SIZE(svfsurf,2)) )
+             tmp(:,:) = svfsurf(:,:)
+             CALL wrd_mpi_io_surface( 'svfsurf', tmp )
+             DEALLOCATE( tmp )
+
           ENDIF
-          IF ( nmrtf > 0 )  THEN
-              WRITE( 89 ) mrtf
-              WRITE( 89 ) mrtft
-              WRITE( 89 ) mrtfsurf
+
+       ENDIF
+
+       IF ( plant_canopy  )  THEN
+
+          lo_no = 0
+          IF ( ncsfl > 0 )  THEN
+
+             DO  j = 1, ncsfl
+                ipcgb = csfsurf(1, j)
+                jc = pcbl(iy,ipcgb)
+                ic = pcbl(ix,ipcgb)
+                lo_no(jc,ic) = lo_no(jc,ic) + 1
+             ENDDO
+
+             ind = 1
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   start_index(j,i) = ind
+                   end_index(j,i)   = start_index(j,i) + lo_no(j,i) - 1
+                   ind = ind+lo_no(j,i)
+                ENDDO
+             ENDDO
+
+          ELSE
+
+             start_index = 1
+             end_index = -1
+
           ENDIF
-!
-!--       Close binary file
-          CALL close_file( 89 )
-
-       ENDIF
-#if defined( __parallel )
-       CALL MPI_BARRIER( comm2d, ierr )
-#endif
-    ENDDO
-
-    CALL location_message( 'writing view factors for radiation interaction', 'finished' )
-
+
+          IF ( ncsfl_tot > 0 )  THEN
+
+             CALL rd_mpi_io_surface_filetypes( start_index, end_index, data_to_write, global_start,&
+                                               global_end )
+             CALL wrd_mpi_io( 'ncsfl_global_start', global_start )
+             CALL wrd_mpi_io( 'ncsfl_global_end', global_end )
+
+             IF ( data_to_write )  THEN
+
+                IF ( ALLOCATED( csf ) )  THEN
+                   CALL wrd_mpi_io_surface( 'csf', csf )
+                ELSE
+                   ALLOCATE( tmp(ndcsf,0) )
+                   CALL wrd_mpi_io_surface( 'csf', tmp )
+                   DEALLOCATE( tmp )
+                ENDIF
+
+                IF ( ALLOCATED( csfsurf ) )  THEN
+                   ALLOCATE( tmp(SIZE(csfsurf,1),SIZE(csfsurf,2)) )
+                   tmp(:,:) = csfsurf(:,:)
+                   CALL wrd_mpi_io_surface( 'csfsurf', tmp )
+                   DEALLOCATE( tmp )
+                ELSE
+                   ALLOCATE( tmp(idcsf,ncsfl) )
+                   CALL wrd_mpi_io_surface( 'csfsurf', tmp )
+                   DEALLOCATE( tmp )
+                ENDIF
+
+             ENDIF
+
+          ENDIF
+
+       ENDIF
+
+       IF ( nmrtbl > 0 )  THEN
+
+          lo_no = 0
+          DO  j = 1, nmrtbl
+             jc = mrtbl(iy,j)
+             ic = mrtbl(ix,j)
+             lo_no(jc,ic) = lo_no(jc,ic) + 1
+          ENDDO
+
+          ind = 1
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                start_index(j,i) = ind
+                end_index(j,i)   = start_index(j,i) + lo_no(j,i) - 1
+                ind = ind+lo_no(j,i)
+             ENDDO
+          ENDDO
+
+       ELSE
+
+          start_index = 1
+          end_index = 0
+       ENDIF
+
+       IF ( nmrtbl_tot > 0 )  THEN
+
+          CALL rd_mpi_io_surface_filetypes( start_index, end_index, data_to_write, global_start,   &
+                                            global_end )
+          CALL wrd_mpi_io( 'nmrtbl_global_start', global_start )
+          CALL wrd_mpi_io( 'nmrtbl_global_end', global_end )
+
+          IF ( data_to_write )  THEN
+
+             CALL wrd_mpi_io_surface( 'mrtsky', mrtsky )
+             CALL wrd_mpi_io_surface( 'mrtskyt', mrtskyt )
+             ALLOCATE( tmp(SIZE(mrtdsit,2),SIZE(mrtdsit,1)) )
+             tmp = TRANSPOSE( mrtdsit )
+             CALL wrd_mpi_io_surface( 'mrtdsit', tmp )
+             DEALLOCATE( tmp )
+
+          ENDIF
+
+       ENDIF
+
+       IF ( nmrtf > 0 )  THEN
+
+          lo_no = 0
+          DO  j = 1, nmrtf
+             isurf = mrtfsurf(1,j)
+             jc = surfl(iy,isurf)
+             ic = surfl(ix,isurf)
+             lo_no(jc,ic) = lo_no(jc,ic) + 1
+          ENDDO
+
+          ind = 1
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                start_index(j,i) = ind
+                end_index(j,i)   = start_index(j,i) + lo_no(j,i) - 1
+                ind = ind+lo_no(j,i)
+             ENDDO
+          ENDDO
+
+       ELSE
+          start_index = 1
+          end_index = 0
+       ENDIF
+
+       IF ( nmrtf_tot > 0 )  THEN
+
+          CALL rd_mpi_io_surface_filetypes( start_index, end_index, data_to_write, global_start,   &
+                                            global_end)
+          CALL wrd_mpi_io( 'nmrtf_global_start', global_start )
+          CALL wrd_mpi_io( 'nmrtf_global_end', global_end )
+
+          IF ( data_to_write )  THEN
+
+             CALL wrd_mpi_io_surface ( 'mrtf', mrtf )
+             CALL wrd_mpi_io_surface ( 'mrtft', mrtft )
+
+             ALLOCATE( tmp(SIZE(mrtfsurf,1),SIZE(mrtfsurf,2)) )
+             tmp = mrtfsurf
+             CALL wrd_mpi_io_surface( 'mrtfsurf', tmp )
+             DEALLOCATE( tmp )
+
+          ENDIF
+
+       ENDIF
+
+       CALL rd_mpi_io_close
+
+       CALL location_message( 'writing sky view factors for radiation interaction ' //             &
+                              '(MPI-IO format)', 'finished' )
+    ENDIF
 
  END SUBROUTINE radiation_write_svf
