Index: palm/trunk/SCRIPTS/.mrun.config.hlrnIII
===================================================================
--- palm/trunk/SCRIPTS/.mrun.config.hlrnIII	(revision 1307)
+++ palm/trunk/SCRIPTS/.mrun.config.hlrnIII	(revision 1308)
@@ -37,11 +37,11 @@
 %compiler_name       ftn                         lccrayb parallel
 %compiler_name_ser   ftn                         lccrayb parallel
-%cpp_options         -e:Z:-DMPI_REAL=MPI_DOUBLE_PRECISION:-DMPI_2REAL=MPI_2DOUBLE_PRECISION:-D__netcdf:-D__netcdf4:-D__fftw   lccrayb parallel
+%cpp_options         -e:Z:-DMPI_REAL=MPI_DOUBLE_PRECISION:-DMPI_2REAL=MPI_2DOUBLE_PRECISION:-D__netcdf:-D__netcdf4:-D__netcdf4_parallel:-D__fftw   lccrayb parallel
 %mopts               -j:4                        lccrayb parallel
 %fopts               -em:-s:real64:-O3:-hnoomp:-hfp3:-hdynamic          lccrayb parallel
 %lopts               -em:-s:real64:-O3:-hnoomp:-hfp3:-hdynamic:-dynamic lccrayb parallel
-%remote_username     <replace by your HLRN-III username>                    lccrayb parallel
-%memory              1500                        lccrayb parallel
-%modules             fftw:cray-netcdf            lccrayb parallel
+%remote_username     <replace by your HLRN-III username>                lccrayb parallel
+%memory              1500                                               lccrayb parallel
+%modules             fftw:cray-hdf5-parallel:cray-netcdf-hdf5parallel   lccrayb parallel
 #
 # HLRN-III Hannover
@@ -50,11 +50,11 @@
 %compiler_name       ftn                         lccrayh parallel
 %compiler_name_ser   ftn                         lccrayh parallel
-%cpp_options         -e:Z:-DMPI_REAL=MPI_DOUBLE_PRECISION:-DMPI_2REAL=MPI_2DOUBLE_PRECISION:-D__netcdf:-D__netcdf4:-D__fftw   lccrayh parallel
+%cpp_options         -e:Z:-DMPI_REAL=MPI_DOUBLE_PRECISION:-DMPI_2REAL=MPI_2DOUBLE_PRECISION:-D__netcdf:-D__netcdf4:-D__netcdf4_parallel:-D__fftw   lccrayh parallel
 %mopts               -j:4                        lccrayh parallel
 %fopts               -em:-s:real64:-O3:-hnoomp:-hfp3:-hdynamic          lccrayh parallel
 %lopts               -em:-s:real64:-O3:-hnoomp:-hfp3:-hdynamic:-dynamic lccrayh parallel
-%remote_username     <replace by your HLRN-III username>                    lccrayh parallel
-%memory              1500                        lccrayh parallel
-%modules             fftw:cray-netcdf            lccrayh parallel
+%remote_username     <replace by your HLRN-III username>                lccrayh parallel
+%memory              1500                                               lccrayh parallel
+%modules             fftw:cray-hdf5-parallel:cray-netcdf-hdf5parallel   lccrayh parallel
 #             
 %write_binary        true                             restart
@@ -70,8 +70,16 @@
 IC:[[ \$localhost = lccrayb ]]  &&  ulimit -v unlimited
 IC:[[ \$localhost = lccrayb ]]  &&  ulimit -s unlimited
+IC:[[ \$localhost = lccrayb ]]  &&  export MPICH_MPIIO_HINTS_DISPLAY=1
+IC:[[ \$localhost = lccrayb ]]  &&  export MPICH_MPIIO_STATS=1
+IC:[[ \$localhost = lccrayb ]]  &&  E8="=8"
+IC:[[ \$localhost = lccrayb ]]  &&  export MPICH_MPIIO_HINTS="*DATA_**:striping_factor"\$E8
 IC:[[ \$localhost = lccrayh ]]  &&  export APRUN_XFER_LIMITS=1
 IC:[[ \$localhost = lccrayh ]]  &&  export APRUN_XFER_STACK_LIMIT=0
 IC:[[ \$localhost = lccrayh ]]  &&  ulimit -v unlimited
 IC:[[ \$localhost = lccrayh ]]  &&  ulimit -s unlimited
+IC:[[ \$localhost = lccrayh ]]  &&  export MPICH_MPIIO_HINTS_DISPLAY=1
+IC:[[ \$localhost = lccrayh ]]  &&  export MPICH_MPIIO_STATS=1
+IC:[[ \$localhost = lccrayh ]]  &&  E8="=8"
+IC:[[ \$localhost = lccrayh ]]  &&  export MPICH_MPIIO_HINTS="*DATA_**:striping_factor"\$E8
 #
 #----------------------------------------------------------------------------
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 1307)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 1308)
@@ -20,5 +20,8 @@
 ! Current revisions:
 ! -----------------
-!
+! +netcdf_data_format_save
+! Calculate fixed number of output time levels for parallel netcdf output.
+! For masked data, parallel netcdf output is not tested so far, hence
+! netcdf_data_format is switched back to non-paralell output.
 !
 ! Former revisions:
@@ -337,5 +340,6 @@
     CHARACTER (LEN=100) ::  action
 
-    INTEGER ::  i, ilen, iremote = 0, j, k, kk, position, prec
+    INTEGER ::  i, ilen, iremote = 0, j, k, kk, netcdf_data_format_save, &
+                position, prec
     LOGICAL ::  found, ldum
     REAL    ::  gradient, remote = 0.0, simulation_time_since_reference
@@ -3265,5 +3269,20 @@
 !
 !--    Generate masks for masked data output
+!--    Parallel netcdf output is not tested so far for masked data, hence 
+!--    netcdf_data_format is switched back to non-paralell output.
+       netcdf_data_format_save = netcdf_data_format
+       IF ( netcdf_data_format > 4 )  THEN
+          IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
+          IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
+          message_string = 'netCDF file formats '//                            &
+                           '5 (parallel netCDF 4) and ' //                     &
+                           '6 (parallel netCDF 4 Classic model) '//            &
+                           '&are currently not supported (not yet tested) ' // &
+                           'for masked data.&Using respective non-parallel' // &  
+                           ' output for masked data.'
+          CALL message( 'check_parameters', 'PA0383', 0, 0, 0, 6, 0 )
+       ENDIF
        CALL init_masks
+       netcdf_data_format = netcdf_data_format_save
     ENDIF
 
@@ -3294,4 +3313,28 @@
     ENDIF
 #endif
+
+!
+!-- Calculate fixed number of output time levels for parallel netcdf output.
+!-- The time dimension has to be limited for paralell output, otherwise the
+!-- performance drops significantly.
+    IF ( netcdf_output )  THEN
+       IF ( netcdf_data_format > 4 )  THEN
+          ntdim_3d(0) = INT( ( end_time - skip_time_do3d ) / dt_do3d )
+          IF ( do3d_at_begin ) ntdim_3d(0) = ntdim_3d(0) + 1
+          ntdim_3d(1) = INT( ( end_time - skip_time_data_output_av ) &
+                              / dt_data_output_av )
+          ntdim_2d_xy(0) = INT( ( end_time - skip_time_do2d_xy ) / dt_do2d_xy )
+          ntdim_2d_xz(0) = INT( ( end_time - skip_time_do2d_xz ) / dt_do2d_xz )
+          ntdim_2d_yz(0) = INT( ( end_time - skip_time_do2d_yz ) / dt_do2d_yz )
+          IF ( do2d_at_begin )  THEN
+             ntdim_2d_xy(0) = ntdim_2d_xy(0) + 1
+             ntdim_2d_xz(0) = ntdim_2d_xz(0) + 1
+             ntdim_2d_yz(0) = ntdim_2d_yz(0) + 1
+          ENDIF
+          ntdim_2d_xy(1) = ntdim_3d(1)
+          ntdim_2d_xz(1) = ntdim_3d(1)
+          ntdim_2d_yz(1) = ntdim_3d(1)
+       ENDIF
+    ENDIF
 
 #if ! defined( __check )
Index: palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 1307)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 1308)
@@ -20,5 +20,10 @@
 ! Current revisions:
 ! -----------------
-!
+! +local_2d_sections, local_2d_sections_l, ns
+! Check, if the limit of the time dimension is exceeded for parallel output
+! To increase the performance for parallel output, the following is done:
+! - Update of time axis is only done by PE0
+! - Cross sections are first stored on a local array and are written
+!   collectively to the output file by all PEs.
 !
 ! Former revisions:
@@ -143,6 +148,6 @@
     CHARACTER (LEN=25) ::  section_chr
     CHARACTER (LEN=50) ::  rtext
-    INTEGER ::  av, ngp, file_id, i, if, is, iis, j, k, l, layer_xy, n, psi, &
-                s, sender, &
+    INTEGER ::  av, ngp, file_id, i, if, is, iis, j, k, l, layer_xy, n, ns, &
+                psi, s, sender, &
                 ind(4)
     LOGICAL ::  found, resorted, two_d
@@ -150,5 +155,6 @@
     REAL, DIMENSION(:), ALLOCATABLE ::      level_z
     REAL, DIMENSION(:,:), ALLOCATABLE ::    local_2d, local_2d_l
-    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  local_pf
+    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  local_pf, local_2d_sections, &
+                                            local_2d_sections_l
 #if defined( __parallel )
     REAL, DIMENSION(:,:),   ALLOCATABLE ::  total_2d
@@ -157,6 +163,4 @@
 
     NAMELIST /LOCAL/  rtext
-
-    CALL cpu_log (log_point(3),'data_output_2d','start')
 
 !
@@ -166,4 +170,36 @@
     IF ( mode == 'xz'  .AND.  .NOT. data_output_xz(av) )  RETURN
     IF ( mode == 'yz'  .AND.  .NOT. data_output_yz(av) )  RETURN
+!
+!-- For parallel netcdf output the time axis must be limited. Return, if this
+!-- limit is exceeded. This could be the case, if the simulated time exceeds 
+!-- the given end time by the length of the given output interval.
+    IF ( netcdf_data_format > 4 )  THEN
+       IF ( mode == 'xy'  .AND.  do2d_xy_time_count(av) + 1 > &
+            ntdim_2d_xy(av) )  THEN
+          WRITE ( message_string, * ) 'Output of xy cross-sections is not ', &
+                          'given at t=', simulated_time, '&because the', & 
+                          ' maximum number of output time levels is exceeded.'
+          CALL message( 'data_output_2d', 'PA0384', 0, 1, 0, 6, 0 )          
+          RETURN
+       ENDIF
+       IF ( mode == 'xz'  .AND.  do2d_xz_time_count(av) + 1 > &
+            ntdim_2d_xz(av) )  THEN
+          WRITE ( message_string, * ) 'Output of xz cross-sections is not ',  &
+                          'given at t=', simulated_time, '&because the', & 
+                          ' maximum number of output time levels is exceeded.'
+          CALL message( 'data_output_2d', 'PA0385', 0, 1, 0, 6, 0 )          
+          RETURN
+       ENDIF
+       IF ( mode == 'yz'  .AND.  do2d_yz_time_count(av) + 1 > &
+            ntdim_2d_yz(av) )  THEN
+          WRITE ( message_string, * ) 'Output of yz cross-sections is not ', &
+                          'given at t=', simulated_time, '&because the', & 
+                          ' maximum number of output time levels is exceeded.'
+          CALL message( 'data_output_2d', 'PA0386', 0, 1, 0, 6, 0 )          
+          RETURN
+       ENDIF
+    ENDIF
+
+    CALL cpu_log (log_point(3),'data_output_2d','start')
 
     two_d = .FALSE.    ! local variable to distinguish between output of pure 2D
@@ -178,4 +214,14 @@
           s = 1
           ALLOCATE( level_z(nzb:nzt+1), local_2d(nxlg:nxrg,nysg:nyng) )
+
+          IF ( netcdf_data_format > 4 )  THEN
+             ns = 1
+             DO WHILE ( section(ns,s) /= -9999  .AND.  ns <= 100 )
+                ns = ns + 1
+             ENDDO
+             ns = ns - 1
+             ALLOCATE( local_2d_sections(nxlg:nxrg,nysg:nyng,1:ns) )
+             local_2d_sections = 0.0
+          ENDIF
 
 !
@@ -201,4 +247,15 @@
           ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
 
+          IF ( netcdf_data_format > 4 )  THEN
+             ns = 1
+             DO WHILE ( section(ns,s) /= -9999  .AND.  ns <= 100 )
+                ns = ns + 1
+             ENDDO
+             ns = ns - 1
+             ALLOCATE( local_2d_sections(nxlg:nxrg,1:ns,nzb:nzt+1) )
+             ALLOCATE( local_2d_sections_l(nxlg:nxrg,1:ns,nzb:nzt+1) )
+             local_2d_sections = 0.0; local_2d_sections_l = 0.0
+          ENDIF
+
 !
 !--       Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
@@ -220,7 +277,17 @@
 
        CASE ( 'yz' )
-
           s = 3
           ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
+
+          IF ( netcdf_data_format > 4 )  THEN
+             ns = 1
+             DO WHILE ( section(ns,s) /= -9999  .AND.  ns <= 100 )
+                ns = ns + 1
+             ENDDO
+             ns = ns - 1
+             ALLOCATE( local_2d_sections(1:ns,nysg:nyng,nzb:nzt+1) )
+             ALLOCATE( local_2d_sections_l(1:ns,nysg:nyng,nzb:nzt+1) )
+             local_2d_sections = 0.0; local_2d_sections_l = 0.0
+          ENDIF
 
 !
@@ -243,5 +310,4 @@
 
        CASE DEFAULT
-
           message_string = 'unknown cross-section: ' // TRIM( mode )
           CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
@@ -799,9 +865,11 @@
 
 !
-!--                Update the netCDF xy cross section time axis
-                   IF ( myid == 0  .OR.  netcdf_data_format > 4 )  THEN
-                      IF ( simulated_time /= do2d_xy_last_time(av) )  THEN
-                         do2d_xy_time_count(av) = do2d_xy_time_count(av) + 1
-                         do2d_xy_last_time(av)  = simulated_time
+!--                Update the netCDF xy cross section time axis.
+!--                In case of parallel output, this is only done by PE0
+!--                to increase the performance.
+                   IF ( simulated_time /= do2d_xy_last_time(av) )  THEN
+                      do2d_xy_time_count(av) = do2d_xy_time_count(av) + 1
+                      do2d_xy_last_time(av)  = simulated_time
+                      IF ( myid == 0 )  THEN
                          IF ( ( .NOT. data_output_2d_on_each_pe  .AND. &
                               netcdf_output )  .OR.  netcdf_data_format > 4 ) &
@@ -846,6 +914,4 @@
 !
 !--                   Parallel output in netCDF4/HDF5 format.
-!--                   Do not output redundant ghost point data except for the
-!--                   boundaries of the total domain.
                       IF ( two_d ) THEN
                          iis = 1
@@ -855,39 +921,13 @@
 
 #if defined( __netcdf )
-                      IF ( nxr == nx  .AND.  nyn /= ny )  THEN
-                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
-                                                 id_var_do2d(av,if),           &
-                                                 local_2d(nxl:nxr+1,nys:nyn),  &
-                                                 start = (/ nxl+1, nys+1, iis, &
-                                                    do2d_xy_time_count(av) /), &
-                                                 count = (/ nxr-nxl+2,         &
-                                                            nyn-nys+1, 1, 1 /) )
-                      ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
-                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
-                                                 id_var_do2d(av,if),           &
-                                                 local_2d(nxl:nxr,nys:nyn+1),  &
-                                                 start = (/ nxl+1, nys+1, iis, &
-                                                    do2d_xy_time_count(av) /), &
-                                                 count = (/ nxr-nxl+1,         &
-                                                            nyn-nys+2, 1, 1 /) )
-                      ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
-                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
-                                                 id_var_do2d(av,if),           &
-                                                 local_2d(nxl:nxr+1,nys:nyn+1),&
-                                                 start = (/ nxl+1, nys+1, iis, &
-                                                    do2d_xy_time_count(av) /), &
-                                                 count = (/ nxr-nxl+2,          &
-                                                            nyn-nys+2, 1, 1 /) )
-                      ELSE
-                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
-                                                 id_var_do2d(av,if),           &
-                                                 local_2d(nxl:nxr,nys:nyn),    &
-                                                 start = (/ nxl+1, nys+1, iis, &
-                                                    do2d_xy_time_count(av) /), &
-                                                 count = (/ nxr-nxl+1,         &
-                                                            nyn-nys+1, 1, 1 /) )
-                      ENDIF
-
-                      CALL handle_netcdf_error( 'data_output_2d', 55 )
+!
+!--                   For parallel output, all cross sections are first stored
+!--                   here on a local array and will be written to the output
+!--                   file afterwards to increase the performance.
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            local_2d_sections(i,j,iis) = local_2d(i,j)
+                         ENDDO
+                      ENDDO
 #endif
                    ELSE
@@ -1047,10 +1087,11 @@
                 CASE ( 'xz' )
 !
-!--                Update the netCDF xz cross section time axis
-                   IF ( myid == 0  .OR.  netcdf_data_format > 4 )  THEN
-
-                      IF ( simulated_time /= do2d_xz_last_time(av) )  THEN
-                         do2d_xz_time_count(av) = do2d_xz_time_count(av) + 1
-                         do2d_xz_last_time(av)  = simulated_time
+!--                Update the netCDF xz cross section time axis.
+!--                In case of parallel output, this is only done by PE0
+!--                to increase the performance.
+                   IF ( simulated_time /= do2d_xz_last_time(av) )  THEN
+                      do2d_xz_time_count(av) = do2d_xz_time_count(av) + 1
+                      do2d_xz_last_time(av)  = simulated_time
+                      IF ( myid == 0 )  THEN
                          IF ( ( .NOT. data_output_2d_on_each_pe  .AND.        &
                               netcdf_output )  .OR.  netcdf_data_format > 4 ) &
@@ -1066,5 +1107,4 @@
                          ENDIF
                       ENDIF
-
                    ENDIF
 
@@ -1114,35 +1154,4 @@
                    IF ( netcdf_output  .AND.  netcdf_data_format > 4 )  THEN
 !
-!--                   ATTENTION: The following lines are a workaround, because
-!--                              independet output does not work with the
-!--                              current netCDF4 installation. Therefore, data
-!--                              are transferred from PEs having the cross
-!--                              sections to other PEs along y having no cross
-!--                              section data. Some of these data are the
-!--                              output.
-!--                   BEGIN WORKAROUND---------------------------------------
-                      IF ( npey /= 1  .AND.  section(is,s) /= -1)  THEN
-                         ALLOCATE( local_2d_l(nxlg:nxrg,nzb:nzt+1) )
-                         local_2d_l = 0.0
-                         IF ( section(is,s) >= nys .AND. section(is,s) <= nyn )&
-                         THEN
-                            local_2d_l = local_2d
-                         ENDIF
-#if defined( __parallel )
-!
-!--                      Distribute data over all PEs along y
-                         ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
-                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
-                         CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb),            &
-                                             local_2d(nxlg,nzb), ngp,         &
-                                             MPI_REAL, MPI_SUM, comm1dy, ierr )
-#else
-                         local_2d = local_2d_l
-#endif
-                         DEALLOCATE( local_2d_l )
-                      ENDIF
-!--                   END WORKAROUND-----------------------------------------
-
-!
 !--                   Output in netCDF4/HDF5 format.
 !--                   Output only on those PEs where the respective cross
@@ -1152,41 +1161,14 @@
                              section(is,s) <= nyn )  .OR.  &
                            ( section(is,s) == -1  .AND.  myidy == 0 ) )  THEN
-!
-!--                      Do not output redundant ghost point data except for the
-!--                      boundaries of the total domain.
 #if defined( __netcdf )
-                         IF ( nxr == nx )  THEN
-                            nc_stat = NF90_PUT_VAR( id_set_xz(av),             &
-                                                id_var_do2d(av,if),            &
-                                                local_2d(nxl:nxr+1,nzb:nzt+1), &
-                                                start = (/ nxl+1, is, 1,       &
-                                                    do2d_xz_time_count(av) /), &
-                                                count = (/ nxr-nxl+2, 1,       &
-                                                           nzt+2, 1 /) )
-                         ELSE
-                            nc_stat = NF90_PUT_VAR( id_set_xz(av),             &
-                                                id_var_do2d(av,if),            &
-                                                local_2d(nxl:nxr,nzb:nzt+1),   &
-                                                start = (/ nxl+1, is, 1,       &
-                                                    do2d_xz_time_count(av) /), &
-                                                count = (/ nxr-nxl+1, 1,       &
-                                                           nzt+2, 1 /) )
-                         ENDIF
-
-                         CALL handle_netcdf_error( 'data_output_2d', 57 )
-
-                      ELSE
-!
-!--                      Output on other PEs. Only one point is output!!
-!--                      ATTENTION: This is a workaround (see above)!!
-                         IF ( npey /= 1 )  THEN
-                            nc_stat = NF90_PUT_VAR( id_set_xz(av),             &
-                                                    id_var_do2d(av,if),        &
-                                                    local_2d(nxl:nxl,nzb:nzb), &
-                                                    start = (/ nxl+1, is, 1,   &
-                                                    do2d_xz_time_count(av) /), &
-                                                    count = (/ 1, 1, 1, 1 /) )
-                            CALL handle_netcdf_error( 'data_output_2d', 451 )
-                         ENDIF
+!
+!--                      For parallel output, all cross sections are first
+!--                      stored here on a local array and will be written to the
+!--                      output file afterwards to increase the performance.
+                         DO  i = nxlg, nxrg
+                            DO  k = nzb, nzt+1
+                               local_2d_sections_l(i,is,k) = local_2d(i,k)
+                            ENDDO
+                         ENDDO
 #endif
                       ENDIF
@@ -1357,10 +1339,11 @@
                 CASE ( 'yz' )
 !
-!--                Update the netCDF yz cross section time axis
-                   IF ( myid == 0  .OR.  netcdf_data_format > 4 )  THEN
-
-                      IF ( simulated_time /= do2d_yz_last_time(av) )  THEN
-                         do2d_yz_time_count(av) = do2d_yz_time_count(av) + 1
-                         do2d_yz_last_time(av)  = simulated_time
+!--                Update the netCDF yz cross section time axis.
+!--                In case of parallel output, this is only done by PE0
+!--                to increase the performance.
+                   IF ( simulated_time /= do2d_yz_last_time(av) )  THEN
+                      do2d_yz_time_count(av) = do2d_yz_time_count(av) + 1
+                      do2d_yz_last_time(av)  = simulated_time
+                      IF ( myid == 0 )  THEN
                          IF ( ( .NOT. data_output_2d_on_each_pe  .AND.        &
                               netcdf_output )  .OR.  netcdf_data_format > 4 ) &
@@ -1376,6 +1359,6 @@
                          ENDIF
                       ENDIF
-
-                   ENDIF
+                   ENDIF
+
 !
 !--                If required, carry out averaging along x
@@ -1423,35 +1406,4 @@
                    IF ( netcdf_output  .AND.  netcdf_data_format > 4 )  THEN
 !
-!--                   ATTENTION: The following lines are a workaround, because
-!--                              independet output does not work with the
-!--                              current netCDF4 installation. Therefore, data
-!--                              are transferred from PEs having the cross
-!--                              sections to other PEs along y having no cross
-!--                              section data. Some of these data are the
-!--                              output.
-!--                   BEGIN WORKAROUND---------------------------------------
-                      IF ( npex /= 1  .AND.  section(is,s) /= -1)  THEN
-                         ALLOCATE( local_2d_l(nysg:nyng,nzb:nzt+1) )
-                         local_2d_l = 0.0
-                         IF ( section(is,s) >= nxl .AND. section(is,s) <= nxr )&
-                         THEN
-                            local_2d_l = local_2d
-                         ENDIF
-#if defined( __parallel )
-!
-!--                      Distribute data over all PEs along x
-                         ngp = ( nyng-nysg+1 ) * ( nzt-nzb + 2 )
-                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
-                         CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb),            &
-                                             local_2d(nysg,nzb), ngp,         &
-                                             MPI_REAL, MPI_SUM, comm1dx, ierr )
-#else
-                         local_2d = local_2d_l
-#endif
-                         DEALLOCATE( local_2d_l )
-                      ENDIF
-!--                   END WORKAROUND-----------------------------------------
-
-!
 !--                   Output in netCDF4/HDF5 format.
 !--                   Output only on those PEs where the respective cross
@@ -1461,41 +1413,14 @@
                              section(is,s) <= nxr )  .OR.  &
                            ( section(is,s) == -1  .AND.  myidx == 0 ) )  THEN
-!
-!--                      Do not output redundant ghost point data except for the
-!--                      boundaries of the total domain.
 #if defined( __netcdf )
-                         IF ( nyn == ny )  THEN
-                            nc_stat = NF90_PUT_VAR( id_set_yz(av),             &
-                                                id_var_do2d(av,if),            &
-                                                local_2d(nys:nyn+1,nzb:nzt+1), &
-                                                start = (/ is, nys+1, 1,       &
-                                                    do2d_yz_time_count(av) /), &
-                                                count = (/ 1, nyn-nys+2,       &
-                                                           nzt+2, 1 /) )
-                         ELSE
-                            nc_stat = NF90_PUT_VAR( id_set_yz(av),             &
-                                                id_var_do2d(av,if),            &
-                                                local_2d(nys:nyn,nzb:nzt+1),   &
-                                                start = (/ is, nys+1, 1,       &
-                                                    do2d_yz_time_count(av) /), &
-                                                count = (/ 1, nyn-nys+1,       &
-                                                           nzt+2, 1 /) )
-                         ENDIF
-
-                         CALL handle_netcdf_error( 'data_output_2d', 60 )
-
-                      ELSE
-!
-!--                      Output on other PEs. Only one point is output!!
-!--                      ATTENTION: This is a workaround (see above)!!
-                         IF ( npex /= 1 )  THEN
-                            nc_stat = NF90_PUT_VAR( id_set_yz(av),             &
-                                                    id_var_do2d(av,if),        &
-                                                    local_2d(nys:nys,nzb:nzb), &
-                                                    start = (/ is, nys+1, 1,   &
-                                                    do2d_yz_time_count(av) /), &
-                                                    count = (/ 1, 1, 1, 1 /) )
-                            CALL handle_netcdf_error( 'data_output_2d', 452 )
-                         ENDIF
+!
+!--                      For parallel output, all cross sections are first
+!--                      stored here on a local array and will be written to the
+!--                      output file afterwards to increase the performance.
+                         DO  j = nysg, nyng
+                            DO  k = nzb, nzt+1
+                               local_2d_sections_l(is,j,k) = local_2d(j,k)
+                            ENDDO
+                         ENDDO
 #endif
                       ENDIF
@@ -1669,4 +1594,174 @@
           ENDDO loop1
 
+!
+!--       For parallel output, all data were collected before on a local array
+!--       and are written now to the netcdf file. This must be done to increase
+!--       the performance of the parallel output.
+#if defined( __netcdf )
+          IF ( netcdf_output .AND. netcdf_data_format > 4 )  THEN
+
+                SELECT CASE ( mode )
+
+                   CASE ( 'xy' )
+                      IF ( two_d ) THEN
+                         iis = 1
+                      ELSE
+                         iis = is-1
+                      ENDIF
+!
+!--                   Do not output redundant ghost point data except for the
+!--                   boundaries of the total domain.
+                      IF ( nxr == nx  .AND.  nyn /= ny )  THEN
+                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
+                                                 id_var_do2d(av,if),           &
+                                                 local_2d_sections(nxl:nxr+1,  &
+                                                    nys:nyn,1:ns),             &
+                                                 start = (/ nxl+1, nys+1, 1,   &
+                                                    do2d_xy_time_count(av) /), &
+                                                 count = (/ nxr-nxl+2,         &
+                                                            nyn-nys+1, ns, 1   &
+                                                          /) )
+                      ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
+                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
+                                                 id_var_do2d(av,if),           &
+                                                 local_2d_sections(nxl:nxr,    &
+                                                    nys:nyn+1,1:ns),           &
+                                                 start = (/ nxl+1, nys+1, 1,   &
+                                                    do2d_xy_time_count(av) /), &
+                                                 count = (/ nxr-nxl+1,         &
+                                                            nyn-nys+2, ns, 1   &
+                                                          /) )
+                      ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
+                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
+                                                 id_var_do2d(av,if),           &
+                                                 local_2d_sections(nxl:nxr+1,  &
+                                                    nys:nyn+1,1:ns),           &
+                                                 start = (/ nxl+1, nys+1, 1,   &
+                                                    do2d_xy_time_count(av) /), &
+                                                 count = (/ nxr-nxl+2,         &
+                                                            nyn-nys+2, ns, 1   &
+                                                          /) )
+                      ELSE
+                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
+                                                 id_var_do2d(av,if),           &
+                                                 local_2d_sections(nxl:nxr,    &
+                                                    nys:nyn,1:ns),             &
+                                                 start = (/ nxl+1, nys+1, 1,   &
+                                                    do2d_xy_time_count(av) /), &
+                                                 count = (/ nxr-nxl+1,         &
+                                                            nyn-nys+1, ns, 1   &
+                                                          /) )
+                      ENDIF    
+
+                      CALL handle_netcdf_error( 'data_output_2d', 55 ) 
+
+                   CASE ( 'xz' )
+!
+!--                   First, all PEs get the information of all cross-sections.
+!--                   Then the data are written to the output file by all PEs
+!--                   while NF90_COLLECTIVE is set in subroutine 
+!--                   define_netcdf_header. Although redundant information are
+!--                   written to the output file in that case, the performance 
+!--                   is significantly better compared to the case where only
+!--                   the first row of PEs in x-direction (myidx = 0) is given
+!--                   the output while NF90_INDEPENDENT is set.
+                      IF ( npey /= 1 )  THEN
+                         
+#if defined( __parallel )
+!
+!--                      Distribute data over all PEs along y
+                         ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 ) * ns
+                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
+                         CALL MPI_ALLREDUCE( local_2d_sections_l(nxlg,1,nzb),  &
+                                             local_2d_sections(nxlg,1,nzb),    &
+                                             ngp, MPI_REAL, MPI_SUM, comm1dy,  &
+                                             ierr )
+#else
+                         local_2d_sections = local_2d_sections_l
+#endif
+                      ENDIF
+!
+!--                   Do not output redundant ghost point data except for the
+!--                   boundaries of the total domain.
+                      IF ( nxr == nx )  THEN
+                         nc_stat = NF90_PUT_VAR( id_set_xz(av),                &
+                                             id_var_do2d(av,if),               & 
+                                             local_2d_sections(nxl:nxr+1,1:ns, &
+                                                nzb:nzt+1),                    &
+                                             start = (/ nxl+1, 1, 1,           &
+                                                do2d_xz_time_count(av) /),     &
+                                             count = (/ nxr-nxl+2, ns, nzt+2,  &
+                                                        1 /) )
+                      ELSE
+                         nc_stat = NF90_PUT_VAR( id_set_xz(av),                &
+                                             id_var_do2d(av,if),               &
+                                             local_2d_sections(nxl:nxr,1:ns,   &
+                                                nzb:nzt+1),                    &
+                                             start = (/ nxl+1, 1, 1,           &
+                                                do2d_xz_time_count(av) /),     &
+                                             count = (/ nxr-nxl+1, ns, nzt+2,  &
+                                                1 /) )
+                      ENDIF
+
+                      CALL handle_netcdf_error( 'data_output_2d', 57 )
+
+                   CASE ( 'yz' )
+!
+!--                   First, all PEs get the information of all cross-sections.
+!--                   Then the data are written to the output file by all PEs
+!--                   while NF90_COLLECTIVE is set in subroutine 
+!--                   define_netcdf_header. Although redundant information are
+!--                   written to the output file in that case, the performance 
+!--                   is significantly better compared to the case where only
+!--                   the first row of PEs in y-direction (myidy = 0) is given
+!--                   the output while NF90_INDEPENDENT is set.
+                      IF ( npex /= 1 )  THEN
+
+#if defined( __parallel )
+!
+!--                      Distribute data over all PEs along x
+                         ngp = ( nyng-nysg+1 ) * ( nzt-nzb + 2 ) * ns
+                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
+                         CALL MPI_ALLREDUCE( local_2d_sections_l(1,nysg,nzb),  &
+                                             local_2d_sections(1,nysg,nzb),    &
+                                             ngp, MPI_REAL, MPI_SUM, comm1dx,  &
+                                             ierr )
+#else
+                         local_2d_sections = local_2d_sections_l
+#endif
+                      ENDIF
+!
+!--                   Do not output redundant ghost point data except for the
+!--                   boundaries of the total domain.
+                      IF ( nyn == ny )  THEN
+                         nc_stat = NF90_PUT_VAR( id_set_yz(av),                &
+                                             id_var_do2d(av,if),               &
+                                             local_2d_sections(1:ns,           &
+                                                nys:nyn+1,nzb:nzt+1),          &
+                                             start = (/ 1, nys+1, 1,           &
+                                                do2d_yz_time_count(av) /),     &
+                                             count = (/ ns, nyn-nys+2,         &
+                                                        nzt+2, 1 /) )
+                      ELSE
+                         nc_stat = NF90_PUT_VAR( id_set_yz(av),                &
+                                             id_var_do2d(av,if),               &
+                                             local_2d_sections(1:ns,nys:nyn,   &
+                                                nzb:nzt+1),                    &
+                                             start = (/ 1, nys+1, 1,           &
+                                                do2d_yz_time_count(av) /),     &
+                                             count = (/ ns, nyn-nys+1,         &
+                                                        nzt+2, 1 /) )
+                      ENDIF
+
+                      CALL handle_netcdf_error( 'data_output_2d', 60 )
+
+                   CASE DEFAULT
+                      message_string = 'unknown cross-section: ' // TRIM( mode )
+                      CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
+
+                END SELECT                     
+
+          ENDIF
+
        ENDIF
 
@@ -1680,5 +1775,8 @@
 !-- Deallocate temporary arrays.
     IF ( ALLOCATED( level_z ) )  DEALLOCATE( level_z )
-    DEALLOCATE( local_pf, local_2d )
+    IF ( netcdf_data_format > 4 )  THEN
+       DEALLOCATE( local_pf, local_2d, local_2d_sections )
+       IF( mode == 'xz' .OR. mode == 'yz' ) DEALLOCATE( local_2d_sections_l )
+    ENDIF
 #if defined( __parallel )
     IF ( .NOT.  data_output_2d_on_each_pe  .AND.  myid == 0 )  THEN
@@ -1704,5 +1802,4 @@
     ENDIF
 
-
     CALL cpu_log (log_point(3),'data_output_2d','stop','nobarrier')
 
Index: palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_3d.f90	(revision 1307)
+++ palm/trunk/SOURCE/data_output_3d.f90	(revision 1308)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! ------------------
-!
+! Check, if the limit of the time dimension is exceeded for parallel output
+! To increase the performance for parallel output, the following is done:
+! - Update of time axis is only done by PE0
 !
 ! Former revisions:
@@ -144,4 +146,18 @@
 !-- Return, if nothing to output
     IF ( do3d_no(av) == 0 )  RETURN
+!
+!-- For parallel netcdf output the time axis must be limited. Return, if this
+!-- limit is exceeded. This could be the case, if the simulated time exceeds 
+!-- the given end time by the length of the given output interval.
+    IF ( netcdf_data_format > 4 )  THEN
+       IF ( do3d_time_count(av) + 1 > ntdim_3d(av) )  THEN
+          WRITE ( message_string, * ) 'Output of 3d data is not given at t=',  &
+                                      simulated_time, '&because the maximum ', & 
+                                      'number of output time levels is ',      &
+                                      'exceeded.'
+          CALL message( 'data_output_3d', 'PA0387', 0, 1, 0, 6, 0 )          
+          RETURN
+       ENDIF
+    ENDIF
 
     CALL cpu_log (log_point(14),'data_output_3d','start')
@@ -174,7 +190,9 @@
 !
 !-- Update the netCDF time axis
+!-- In case of parallel output, this is only done by PE0 to increase the
+!-- performance.
 #if defined( __netcdf )
-    IF ( myid == 0  .OR.  netcdf_data_format > 4 )  THEN
-       do3d_time_count(av) = do3d_time_count(av) + 1
+    do3d_time_count(av) = do3d_time_count(av) + 1
+    IF ( myid == 0 )  THEN
        IF ( netcdf_output )  THEN
           nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_time_3d(av), &
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 1307)
+++ palm/trunk/SOURCE/header.f90	(revision 1308)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-!
+! output of the fixed number of output time levels
+! output_format adjusted for masked data if netcdf_data_format > 5
 !
 ! Former revisions:
@@ -924,5 +925,9 @@
                 ENDIF
              ENDIF
-
+             IF ( netcdf_data_format > 4 )  THEN
+                WRITE ( io, 352 )  ntdim_2d_xy(av)
+             ELSE
+                WRITE ( io, 353 )
+             ENDIF
           ENDIF
 
@@ -962,4 +967,9 @@
                 ENDIF
              ENDIF
+             IF ( netcdf_data_format > 4 )  THEN
+                WRITE ( io, 352 )  ntdim_2d_xz(av)
+             ELSE
+                WRITE ( io, 353 )
+             ENDIF
           ENDIF
 
@@ -999,4 +1009,9 @@
                 ENDIF
              ENDIF
+             IF ( netcdf_data_format > 4 )  THEN
+                WRITE ( io, 352 )  ntdim_2d_yz(av)
+             ELSE
+                WRITE ( io, 353 )
+             ENDIF
           ENDIF
 
@@ -1050,4 +1065,10 @@
                                 TRIM( begin_chr ), averaging_interval, &
                                 dt_averaging_input, zu(nz_do3d), nz_do3d
+          ENDIF
+
+          IF ( netcdf_data_format > 4 )  THEN
+             WRITE ( io, 352 )  ntdim_3d(av)
+          ELSE
+             WRITE ( io, 353 )
           ENDIF
 
@@ -1119,4 +1140,8 @@
                 output_format = output_format_netcdf
              ENDIF
+!--          Parallel output not implemented for mask data, hence
+!--          output_format must be adjusted.
+             IF ( netcdf_data_format == 5 ) output_format = 'netCDF4/HDF5'
+             IF ( netcdf_data_format == 6 ) output_format = 'netCDF4/HDF5 classic'
              WRITE ( io, 344 )  output_format
 
@@ -1908,4 +1933,6 @@
             'mask_scale_',A,' constructed from array mask_',I2.2,'_',A,'_loop:'/ &
             '          loop begin:',F8.2,', end:',F8.2,', stride:',F8.2 )
+352 FORMAT  (/'       Fixed number of output time levels allowed: ',I3 /)
+353 FORMAT  (/'       Fixed number of output time levels allowed: unlimited' /)
 #if defined( __dvrp_graphics )
 360 FORMAT ('    Plot-Sequence with dvrp-software:'/ &
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 1307)
+++ palm/trunk/SOURCE/modules.f90	(revision 1308)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! +ntdim_2d_xy, ntdim_2d_xz, ntdim_2d_yz, ntdim_3d
 !
 ! Former revisions:
@@ -760,4 +760,6 @@
     INTEGER :: check_restart = 0
 #endif
+
+    INTEGER, DIMENSION(0:1) :: ntdim_2d_xy, ntdim_2d_xz, ntdim_2d_yz, ntdim_3d
 
     INTEGER, DIMENSION(:), ALLOCATABLE ::  grid_level_count
Index: palm/trunk/SOURCE/netcdf.f90
===================================================================
--- palm/trunk/SOURCE/netcdf.f90	(revision 1307)
+++ palm/trunk/SOURCE/netcdf.f90	(revision 1308)
@@ -23,5 +23,11 @@
 ! Current revisions:
 ! ------------------
-!
+! +ntime_count, oldmode
+! Adjust NF90_CREATE and NF90_OPEN statement for parallel output
+! To increase the performance for parallel output, the following is done:
+! - Limit time dimension
+! - Values of axis data are only written by PE0
+! - No fill is set for all variables
+! Check the number of output time levels for restart jobs
 !
 ! Former revisions:
@@ -169,5 +175,7 @@
     INTEGER ::  av, cross_profiles_count, cross_profiles_maxi, delim, &
                 delim_old, file_id, i, id_last, id_x, id_y, id_z, j,  &
-                k, kk, ns, ns_old, nz_old
+                k, kk, ns, ns_old, ntime_count, nz_old
+
+    INTEGER, SAVE ::  oldmode
 
     INTEGER, DIMENSION(1) ::  id_dim_time_old, id_dim_x_yz_old,  &
@@ -858,8 +866,16 @@
 
 !
-!--       Define time coordinate for volume data (unlimited dimension)
-          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'time', NF90_UNLIMITED, &
-                                  id_dim_time_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 64 )
+!--       Define time coordinate for volume data.
+!--       For parallel output the time dimensions has to be limited, otherwise
+!--       the performance drops significantly.
+          IF ( netcdf_data_format < 5 )  THEN
+             nc_stat = NF90_DEF_DIM( id_set_3d(av), 'time', NF90_UNLIMITED, &
+                                     id_dim_time_3d(av) )
+             CALL handle_netcdf_error( 'netcdf', 64 )
+          ELSE
+             nc_stat = NF90_DEF_DIM( id_set_3d(av), 'time', ntdim_3d(av), &
+                                     id_dim_time_3d(av) )
+             CALL handle_netcdf_error( 'netcdf', 523 )
+          ENDIF
 
           nc_stat = NF90_DEF_VAR( id_set_3d(av), 'time', NF90_DOUBLE, &
@@ -1088,7 +1104,13 @@
              CALL handle_netcdf_error( 'netcdf', 357 )
 #if defined( __netcdf4_parallel )
-!
-!--          Set collective io operations for parallel io
              IF ( netcdf_data_format > 4 )  THEN
+!
+!--             Set no fill for every variable to increase performance.
+                nc_stat = NF90_DEF_VAR_FILL( id_set_3d(av),     &
+                                             id_var_do3d(av,i), &
+                                             1, 0 )
+                CALL handle_netcdf_error( 'netcdf', 532 )
+!
+!--             Set collective io operations for parallel io
                 nc_stat = NF90_VAR_PAR_ACCESS( id_set_3d(av),     &
                                                id_var_do3d(av,i), &
@@ -1113,4 +1135,10 @@
 
 !
+!--       Set general no fill, otherwise the performance drops significantly for
+!--       parallel output.
+          nc_stat = NF90_SET_FILL( id_set_3d(av), NF90_NOFILL, oldmode )
+          CALL handle_netcdf_error( 'netcdf', 528 )
+
+!
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_3d(av) )
@@ -1118,79 +1146,85 @@
 
 !
-!--       Write data for x (shifted by +dx/2) and xu axis
-          ALLOCATE( netcdf_data(0:nx+1) )
-
-          DO  i = 0, nx+1
-             netcdf_data(i) = ( i + 0.5 ) * dx
-          ENDDO
-
-          nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_x_3d(av), netcdf_data, &
-                                  start = (/ 1 /), count = (/ nx+2 /) )
-          CALL handle_netcdf_error( 'netcdf', 83 )
-
-          DO  i = 0, nx+1
-             netcdf_data(i) = i * dx
-          ENDDO
-
-          nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_xu_3d(av), &
-                                  netcdf_data, start = (/ 1 /),    &
-                                  count = (/ nx+2 /) )
-          CALL handle_netcdf_error( 'netcdf', 385 )
-
-          DEALLOCATE( netcdf_data )
-
-!
-!--       Write data for y (shifted by +dy/2) and yv axis
-          ALLOCATE( netcdf_data(0:ny+1) )
-
-          DO  i = 0, ny+1
-             netcdf_data(i) = ( i + 0.5 ) * dy
-          ENDDO
-
-          nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_y_3d(av), netcdf_data, &
-                                  start = (/ 1 /), count = (/ ny+2 /))
-          CALL handle_netcdf_error( 'netcdf', 84 )
-
-          DO  i = 0, ny+1
-             netcdf_data(i) = i * dy
-          ENDDO
-
-          nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_yv_3d(av), &
-                                  netcdf_data, start = (/ 1 /),    &
-                                  count = (/ ny+2 /))
-          CALL handle_netcdf_error( 'netcdf', 387 )
-
-          DEALLOCATE( netcdf_data )
-
-!
-!--       Write zu and zw data (vertical axes)
-          nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zu_3d(av),    &
-                                  zu(nzb:nz_do3d), start = (/ 1 /), &
-                                  count = (/ nz_do3d-nzb+1 /) )
-          CALL handle_netcdf_error( 'netcdf', 85 )
-
-
-          nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zw_3d(av),    &
-                                  zw(nzb:nz_do3d), start = (/ 1 /), &
-                                  count = (/ nz_do3d-nzb+1 /) )
-          CALL handle_netcdf_error( 'netcdf', 86 )
-
-
-!
-!--       In case of non-flat topography write height information
-          IF ( TRIM( topography ) /= 'flat' )  THEN
-
-             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av), &
-                                     zu_s_inner(0:nx+1,0:ny+1), &
-                                     start = (/ 1, 1 /), &
-                                     count = (/ nx+2, ny+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 419 )
-
-             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av), &
-                                     zw_w_inner(0:nx+1,0:ny+1), &
-                                     start = (/ 1, 1 /), &
-                                     count = (/ nx+2, ny+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 420 )
-
+!--       These data are only written by PE0 for parallel output to increase
+!--       the performance.
+          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
+!
+!--          Write data for x (shifted by +dx/2) and xu axis
+             ALLOCATE( netcdf_data(0:nx+1) )
+
+             DO  i = 0, nx+1
+                netcdf_data(i) = ( i + 0.5 ) * dx
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_x_3d(av),  &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ nx+2 /) )
+             CALL handle_netcdf_error( 'netcdf', 83 )
+
+             DO  i = 0, nx+1
+                netcdf_data(i) = i * dx
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_xu_3d(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ nx+2 /) )
+             CALL handle_netcdf_error( 'netcdf', 385 )
+
+             DEALLOCATE( netcdf_data )
+
+!
+!--          Write data for y (shifted by +dy/2) and yv axis
+             ALLOCATE( netcdf_data(0:ny+1) )
+
+             DO  i = 0, ny+1
+                netcdf_data(i) = ( i + 0.5 ) * dy
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_y_3d(av),  &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ ny+2 /) )
+             CALL handle_netcdf_error( 'netcdf', 84 )
+
+             DO  i = 0, ny+1
+                netcdf_data(i) = i * dy
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_yv_3d(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ ny+2 /))
+             CALL handle_netcdf_error( 'netcdf', 387 )
+
+             DEALLOCATE( netcdf_data )
+
+!
+!--          Write zu and zw data (vertical axes)
+             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zu_3d(av),  &
+                                     zu(nzb:nz_do3d), start = (/ 1 /), &
+                                     count = (/ nz_do3d-nzb+1 /) )
+             CALL handle_netcdf_error( 'netcdf', 85 )
+
+
+             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zw_3d(av),  &
+                                     zw(nzb:nz_do3d), start = (/ 1 /), &
+                                     count = (/ nz_do3d-nzb+1 /) )
+             CALL handle_netcdf_error( 'netcdf', 86 )
+
+!
+!--          In case of non-flat topography write height information
+             IF ( TRIM( topography ) /= 'flat' )  THEN
+
+                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av), &
+                                        zu_s_inner(0:nx+1,0:ny+1), &
+                                        start = (/ 1, 1 /), &
+                                        count = (/ nx+2, ny+2 /) )
+                CALL handle_netcdf_error( 'netcdf', 419 )
+
+                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av), &
+                                        zw_w_inner(0:nx+1,0:ny+1), &
+                                        start = (/ 1, 1 /), &
+                                        count = (/ nx+2, ny+2 /) )
+                CALL handle_netcdf_error( 'netcdf', 420 )
+
+             ENDIF
           ENDIF
 
@@ -1271,6 +1305,13 @@
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_time_3d(av), &
-                                            len = do3d_time_count(av) )
+                                            len = ntime_count )
           CALL handle_netcdf_error( 'netcdf', 93 )
+
+!
+!--       For non-parallel output use the last output time level of the netcdf
+!--       file because the time dimension is unlimited. In case of parallel
+!--       output the variable ntime_count could get the value of 9*10E36 because
+!--       the time dimension is limited.
+          IF ( netcdf_data_format < 5 ) do3d_time_count(av) = ntime_count
 
           nc_stat = NF90_GET_VAR( id_set_3d(av), id_var_time_3d(av), &
@@ -1294,4 +1335,20 @@
           ENDIF
 
+          IF ( netcdf_data_format > 4 )  THEN
+             IF ( ntdim_3d(av) > ntime_count )  THEN
+                message_string = 'netCDF file for volume data ' // &
+                                 TRIM( var ) // ' from previous run found,' // &
+                                 '&but this file cannot be extended becaus' // &
+                                 'e the number of output time levels&has b' // &
+                                 'een increased compared to the previous s' // &
+                                 'imulation.' //                                &
+                                 '&New file is created instead.'
+                CALL message( 'define_netcdf_header', 'PA0388', 0, 1, 0, 6, 0 )
+                do3d_time_count(av) = 0
+                extend = .FALSE.
+                RETURN
+             ENDIF
+          ENDIF
+
 !
 !--       Dataset seems to be extendable.
@@ -1367,8 +1424,16 @@
 
 !
-!--       Define time coordinate for xy sections (unlimited dimension)
-          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'time', NF90_UNLIMITED, &
-                                  id_dim_time_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 99 )
+!--       Define time coordinate for xy sections.
+!--       For parallel output the time dimensions has to be limited, otherwise
+!--       the performance drops significantly.
+          IF ( netcdf_data_format < 5 )  THEN
+             nc_stat = NF90_DEF_DIM( id_set_xy(av), 'time', NF90_UNLIMITED, &
+                                     id_dim_time_xy(av) )
+             CALL handle_netcdf_error( 'netcdf', 99 )
+          ELSE
+             nc_stat = NF90_DEF_DIM( id_set_xy(av), 'time', ntdim_2d_xy(av), &
+                                     id_dim_time_xy(av) )
+             CALL handle_netcdf_error( 'netcdf', 524 )
+          ENDIF
 
           nc_stat = NF90_DEF_VAR( id_set_xy(av), 'time', NF90_DOUBLE, &
@@ -1655,7 +1720,13 @@
                 CALL handle_netcdf_error( 'netcdf', 354 ) 
 #if defined( __netcdf4_parallel )
-!
-!--             Set collective io operations for parallel io
                 IF ( netcdf_data_format > 4 )  THEN
+!
+!--                Set no fill for every variable to increase performance.
+                   nc_stat = NF90_DEF_VAR_FILL( id_set_xy(av),     &
+                                                id_var_do2d(av,i), &
+                                                1, 0 )
+                   CALL handle_netcdf_error( 'netcdf', 533 )
+!
+!--                Set collective io operations for parallel io
                    nc_stat = NF90_VAR_PAR_ACCESS( id_set_xy(av),     &
                                                   id_var_do2d(av,i), &
@@ -1682,4 +1753,10 @@
 
 !
+!--       Set general no fill, otherwise the performance drops significantly for
+!--       parallel output.
+          nc_stat = NF90_SET_FILL( id_set_xy(av), NF90_NOFILL, oldmode )
+          CALL handle_netcdf_error( 'netcdf', 529 )
+
+!
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_xy(av) )
@@ -1687,116 +1764,123 @@
 
 !
-!--       Write axis data: z_xy, x, y
-          ALLOCATE( netcdf_data(1:ns) )
-
-!
-!--       Write zu data
-          DO  i = 1, ns
-             IF( section(i,1) == -1 )  THEN
-                netcdf_data(i) = -1.0  ! section averaged along z
-             ELSE
-                netcdf_data(i) = zu( section(i,1) )
+!--       These data are only written by PE0 for parallel output to increase
+!--       the performance.
+          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
+
+!
+!--          Write axis data: z_xy, x, y
+             ALLOCATE( netcdf_data(1:ns) )
+
+!
+!--          Write zu data
+             DO  i = 1, ns
+                IF( section(i,1) == -1 )  THEN
+                   netcdf_data(i) = -1.0  ! section averaged along z
+                ELSE
+                   netcdf_data(i) = zu( section(i,1) )
+                ENDIF
+             ENDDO
+             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu_xy(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ ns /) )
+             CALL handle_netcdf_error( 'netcdf', 123 )
+
+!
+!--          Write zw data
+             DO  i = 1, ns
+                IF( section(i,1) == -1 )  THEN
+                   netcdf_data(i) = -1.0  ! section averaged along z
+                ELSE
+                   netcdf_data(i) = zw( section(i,1) )
+                ENDIF
+             ENDDO
+             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zw_xy(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ ns /) )
+             CALL handle_netcdf_error( 'netcdf', 124 )
+
+!
+!--          Write gridpoint number data
+             netcdf_data(1:ns) = section(1:ns,1)
+             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_ind_z_xy(av), &
+                                     netcdf_data, start = (/ 1 /),       &
+                                     count = (/ ns /) )
+             CALL handle_netcdf_error( 'netcdf', 125 )
+
+             DEALLOCATE( netcdf_data )
+
+!
+!--          Write the cross section height u*, t*
+             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu1_xy(av), &
+                                     (/ zu(nzb+1) /), start = (/ 1 /), &
+                                     count = (/ 1 /) )
+             CALL handle_netcdf_error( 'netcdf', 126 )
+
+!
+!--          Write data for x (shifted by +dx/2) and xu axis
+             ALLOCATE( netcdf_data(0:nx+1) )
+
+             DO  i = 0, nx+1
+                netcdf_data(i) = ( i + 0.5 ) * dx
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_x_xy(av), &
+                                     netcdf_data, start = (/ 1 /),   &
+                                     count = (/ nx+2 /) )
+             CALL handle_netcdf_error( 'netcdf', 127 )
+
+             DO  i = 0, nx+1
+                netcdf_data(i) = i * dx
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_xu_xy(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ nx+2 /) )
+             CALL handle_netcdf_error( 'netcdf', 367 )
+
+             DEALLOCATE( netcdf_data )
+
+!
+!--          Write data for y (shifted by +dy/2) and yv axis
+             ALLOCATE( netcdf_data(0:ny+1) )
+
+             DO  i = 0, ny+1
+                netcdf_data(i) = ( i + 0.5 ) * dy
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_y_xy(av), &
+                                     netcdf_data, start = (/ 1 /),   &
+                                     count = (/ ny+2 /))
+             CALL handle_netcdf_error( 'netcdf', 128 )
+
+             DO  i = 0, ny+1
+                netcdf_data(i) = i * dy
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_yv_xy(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ ny+2 /))
+             CALL handle_netcdf_error( 'netcdf', 368 )
+
+             DEALLOCATE( netcdf_data )
+
+!
+!--          In case of non-flat topography write height information
+             IF ( TRIM( topography ) /= 'flat' )  THEN
+
+                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av), &
+                                        zu_s_inner(0:nx+1,0:ny+1), &
+                                        start = (/ 1, 1 /), &
+                                        count = (/ nx+2, ny+2 /) )
+                CALL handle_netcdf_error( 'netcdf', 427 )
+
+                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av), &
+                                        zw_w_inner(0:nx+1,0:ny+1), &
+                                        start = (/ 1, 1 /), &
+                                        count = (/ nx+2, ny+2 /) )
+                CALL handle_netcdf_error( 'netcdf', 428 )
+
              ENDIF
-          ENDDO
-          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu_xy(av), &
-                                  netcdf_data, start = (/ 1 /),    &
-                                  count = (/ ns /) )
-          CALL handle_netcdf_error( 'netcdf', 123 )
-
-!
-!--       Write zw data
-          DO  i = 1, ns
-             IF( section(i,1) == -1 )  THEN
-                netcdf_data(i) = -1.0  ! section averaged along z
-             ELSE
-                netcdf_data(i) = zw( section(i,1) )
-             ENDIF
-          ENDDO
-          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zw_xy(av), &
-                                  netcdf_data, start = (/ 1 /),    &
-                                  count = (/ ns /) )
-          CALL handle_netcdf_error( 'netcdf', 124 )
-
-!
-!--       Write gridpoint number data
-          netcdf_data(1:ns) = section(1:ns,1)
-          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_ind_z_xy(av), &
-                                  netcdf_data, start = (/ 1 /),       &
-                                  count = (/ ns /) )
-          CALL handle_netcdf_error( 'netcdf', 125 )
-
-          DEALLOCATE( netcdf_data )
-
-!
-!--       Write the cross section height u*, t*
-          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu1_xy(av), &
-                                  (/ zu(nzb+1) /), start = (/ 1 /), &
-                                  count = (/ 1 /) )
-          CALL handle_netcdf_error( 'netcdf', 126 )
-
-!
-!--       Write data for x (shifted by +dx/2) and xu axis
-          ALLOCATE( netcdf_data(0:nx+1) )
-
-          DO  i = 0, nx+1
-             netcdf_data(i) = ( i + 0.5 ) * dx
-          ENDDO
-
-          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_x_xy(av), netcdf_data, &
-                                  start = (/ 1 /), count = (/ nx+2 /) )
-          CALL handle_netcdf_error( 'netcdf', 127 )
-
-          DO  i = 0, nx+1
-             netcdf_data(i) = i * dx
-          ENDDO
-
-          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_xu_xy(av), &
-                                  netcdf_data, start = (/ 1 /),    &
-                                  count = (/ nx+2 /) )
-          CALL handle_netcdf_error( 'netcdf', 367 )
-
-          DEALLOCATE( netcdf_data )
-
-!
-!--       Write data for y (shifted by +dy/2) and yv axis
-          ALLOCATE( netcdf_data(0:ny+1) )
-
-          DO  i = 0, ny+1
-             netcdf_data(i) = ( i + 0.5 ) * dy
-          ENDDO
-
-          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_y_xy(av), netcdf_data, &
-                                  start = (/ 1 /), count = (/ ny+2 /))
-          CALL handle_netcdf_error( 'netcdf', 128 )
-
-          DO  i = 0, ny+1
-             netcdf_data(i) = i * dy
-          ENDDO
-
-          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_yv_xy(av), &
-                                  netcdf_data, start = (/ 1 /),    &
-                                  count = (/ ny+2 /))
-          CALL handle_netcdf_error( 'netcdf', 368 )
-
-          DEALLOCATE( netcdf_data )
-
-!
-!--       In case of non-flat topography write height information
-          IF ( TRIM( topography ) /= 'flat' )  THEN
-
-             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av), &
-                                     zu_s_inner(0:nx+1,0:ny+1), &
-                                     start = (/ 1, 1 /), &
-                                     count = (/ nx+2, ny+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 427 )
-
-             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av), &
-                                     zw_w_inner(0:nx+1,0:ny+1), &
-                                     start = (/ 1, 1 /), &
-                                     count = (/ nx+2, ny+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 428 )
-
-          ENDIF
-
+          ENDIF
 
        CASE ( 'xy_ext' )
@@ -1924,6 +2008,13 @@
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_time_xy(av), &
-                                            len = do2d_xy_time_count(av) )
+                                            len = ntime_count )
           CALL handle_netcdf_error( 'netcdf', 136 )
+
+!
+!--       For non-parallel output use the last output time level of the netcdf
+!--       file because the time dimension is unlimited. In case of parallel
+!--       output the variable ntime_count could get the value of 9*10E36 because
+!--       the time dimension is limited.
+          IF ( netcdf_data_format < 5 ) do2d_xy_time_count(av) = ntime_count
 
           nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_time_xy(av),    &
@@ -1945,4 +2036,20 @@
              extend = .FALSE.
              RETURN
+          ENDIF
+
+          IF ( netcdf_data_format > 4 )  THEN
+             IF ( ntdim_2d_xy(av) > ntime_count )  THEN
+                message_string = 'netCDF file for cross sections ' //          &
+                                 TRIM( var ) // ' from previous run found,' // &
+                                 '&but this file cannot be extended becaus' // &
+                                 'e the number of output time levels&has b' // &
+                                 'een increased compared to the previous s' // &
+                                 'imulation.' //                                &
+                                 '&New file is created instead.'
+                CALL message( 'define_netcdf_header', 'PA0389', 0, 1, 0, 6, 0 )
+                do2d_xy_time_count(av) = 0
+                extend = .FALSE.
+                RETURN
+             ENDIF
           ENDIF
 
@@ -2023,8 +2130,16 @@
 
 !
-!--       Define time coordinate for xz sections (unlimited dimension)
-          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'time', NF90_UNLIMITED, &
-                                  id_dim_time_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 142 )
+!--       Define time coordinate for xz sections.
+!--       For parallel output the time dimensions has to be limited, otherwise
+!--       the performance drops significantly.
+          IF ( netcdf_data_format < 5 )  THEN
+             nc_stat = NF90_DEF_DIM( id_set_xz(av), 'time', NF90_UNLIMITED, &
+                                     id_dim_time_xz(av) )
+             CALL handle_netcdf_error( 'netcdf', 142 )
+          ELSE
+             nc_stat = NF90_DEF_DIM( id_set_xz(av), 'time', ntdim_2d_xz(av), &
+                                     id_dim_time_xz(av) )
+             CALL handle_netcdf_error( 'netcdf', 525 )
+          ENDIF
 
           nc_stat = NF90_DEF_VAR( id_set_xz(av), 'time', NF90_DOUBLE, &
@@ -2240,9 +2355,16 @@
                 CALL handle_netcdf_error( 'netcdf', 355 )
 #if defined( __netcdf4_parallel )
-!
-!--             Set independent io operations for parallel io. Collective io
-!--             is only allowed in case of a 1d-decomposition along x, because
-!--             otherwise, not all PEs have output data.
+
                 IF ( netcdf_data_format > 4 )  THEN
+!
+!--                Set no fill for every variable to increase performance.
+                   nc_stat = NF90_DEF_VAR_FILL( id_set_xz(av),     &
+                                                id_var_do2d(av,i), &
+                                                1, 0 )
+                   CALL handle_netcdf_error( 'netcdf', 534 )
+!
+!--                Set independent io operations for parallel io. Collective io
+!--                is only allowed in case of a 1d-decomposition along x,
+!--                because otherwise, not all PEs have output data.
                    IF ( npey == 1 )  THEN
                       nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
@@ -2251,11 +2373,11 @@
                    ELSE
 !
-!--                   ATTENTION: Due to a probable bug in the netCDF4
-!--                              installation, independet output does not work
-!--                              A workaround is used in data_output_2d on those
-!--                              PEs having no data
-                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
-                                                     id_var_do2d(av,i), &
-                                                     NF90_COLLECTIVE )
+!--                   Test simulations showed that the output of cross sections
+!--                   by all PEs in data_output_2d using NF90_COLLECTIVE is 
+!--                   faster than the output by the first row of PEs in
+!--                   x-direction using NF90_INDEPENDENT.
+                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),    & 
+                                                    id_var_do2d(av,i), &
+                                                    NF90_COLLECTIVE )
 !                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
 !                                                     id_var_do2d(av,i), &
@@ -2283,4 +2405,10 @@
 
 !
+!--       Set general no fill, otherwise the performance drops significantly for
+!--       parallel output.
+          nc_stat = NF90_SET_FILL( id_set_xz(av), NF90_NOFILL, oldmode )
+          CALL handle_netcdf_error( 'netcdf', 530 )
+
+!
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_xz(av) )
@@ -2288,85 +2416,94 @@
 
 !
-!--       Write axis data: y_xz, x, zu, zw
-          ALLOCATE( netcdf_data(1:ns) )
-
-!
-!--       Write y_xz data (shifted by +dy/2)
-          DO  i = 1, ns
-             IF( section(i,2) == -1 )  THEN
-                netcdf_data(i) = -1.0  ! section averaged along y
-             ELSE
-                netcdf_data(i) = ( section(i,2) + 0.5 ) * dy
-             ENDIF
-          ENDDO
-          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_y_xz(av), netcdf_data, &
-                                  start = (/ 1 /), count = (/ ns /) )
-          CALL handle_netcdf_error( 'netcdf', 163 )
-
-!
-!--       Write yv_xz data
-          DO  i = 1, ns
-             IF( section(i,2) == -1 )  THEN
-                netcdf_data(i) = -1.0  ! section averaged along y
-             ELSE
-                netcdf_data(i) = section(i,2) * dy
-             ENDIF
-          ENDDO
-          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_yv_xz(av), &
-                                  netcdf_data, start = (/ 1 /),    &
-                                  count = (/ ns /) )
-          CALL handle_netcdf_error( 'netcdf', 375 )
-
-!
-!--       Write gridpoint number data
-          netcdf_data(1:ns) = section(1:ns,2)
-          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_ind_y_xz(av), &
-                                  netcdf_data, start = (/ 1 /),       &
-                                  count = (/ ns /) )
-          CALL handle_netcdf_error( 'netcdf', 164 )
-
-
-          DEALLOCATE( netcdf_data )
-
-!
-!--       Write data for x (shifted by +dx/2) and xu axis
-          ALLOCATE( netcdf_data(0:nx+1) )
-
-          DO  i = 0, nx+1
-             netcdf_data(i) = ( i + 0.5 ) * dx
-          ENDDO
-
-          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_x_xz(av), netcdf_data, &
-                                  start = (/ 1 /), count = (/ nx+2 /) )
-          CALL handle_netcdf_error( 'netcdf', 165 )
-
-          DO  i = 0, nx+1
-             netcdf_data(i) = i * dx
-          ENDDO
-
-          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_xu_xz(av), &
-                                  netcdf_data, start = (/ 1 /),    &
-                                  count = (/ nx+2 /) )
-          CALL handle_netcdf_error( 'netcdf', 377 )
-
-          DEALLOCATE( netcdf_data )
-
-!
-!--       Write zu and zw data (vertical axes)
-          ALLOCATE( netcdf_data(0:nz+1) )
-
-          netcdf_data(0:nz+1) = zu(nzb:nzt+1)
-          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zu_xz(av), &
-                                  netcdf_data, start = (/ 1 /),    &
-                                  count = (/ nz+2 /) )
-          CALL handle_netcdf_error( 'netcdf', 166 )
-
-          netcdf_data(0:nz+1) = zw(nzb:nzt+1)
-          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zw_xz(av), &
-                                  netcdf_data, start = (/ 1 /),    &
-                                  count = (/ nz+2 /) )
-          CALL handle_netcdf_error( 'netcdf', 167 )
-
-          DEALLOCATE( netcdf_data )
+!--       These data are only written by PE0 for parallel output to increase
+!--       the performance.
+          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
+
+!
+!--          Write axis data: y_xz, x, zu, zw
+             ALLOCATE( netcdf_data(1:ns) )
+
+!
+!--          Write y_xz data (shifted by +dy/2)
+             DO  i = 1, ns
+                IF( section(i,2) == -1 )  THEN
+                   netcdf_data(i) = -1.0  ! section averaged along y
+                ELSE
+                   netcdf_data(i) = ( section(i,2) + 0.5 ) * dy
+                ENDIF
+             ENDDO
+             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_y_xz(av), &
+                                     netcdf_data, start = (/ 1 /),   &
+                                     count = (/ ns /) )
+             CALL handle_netcdf_error( 'netcdf', 163 )
+
+!
+!--          Write yv_xz data
+             DO  i = 1, ns
+                IF( section(i,2) == -1 )  THEN
+                   netcdf_data(i) = -1.0  ! section averaged along y
+                ELSE
+                   netcdf_data(i) = section(i,2) * dy
+                ENDIF
+             ENDDO
+             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_yv_xz(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ ns /) )
+             CALL handle_netcdf_error( 'netcdf', 375 )
+
+!
+!--          Write gridpoint number data
+             netcdf_data(1:ns) = section(1:ns,2)
+             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_ind_y_xz(av), &
+                                     netcdf_data, start = (/ 1 /),       &
+                                     count = (/ ns /) )
+             CALL handle_netcdf_error( 'netcdf', 164 )
+
+
+             DEALLOCATE( netcdf_data )
+
+!
+!--          Write data for x (shifted by +dx/2) and xu axis
+             ALLOCATE( netcdf_data(0:nx+1) )
+
+             DO  i = 0, nx+1
+                netcdf_data(i) = ( i + 0.5 ) * dx
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_x_xz(av), &
+                                     netcdf_data, start = (/ 1 /),   &
+                                     count = (/ nx+2 /) )
+             CALL handle_netcdf_error( 'netcdf', 165 )
+
+             DO  i = 0, nx+1
+                netcdf_data(i) = i * dx
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_xu_xz(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ nx+2 /) )
+             CALL handle_netcdf_error( 'netcdf', 377 )
+
+             DEALLOCATE( netcdf_data )
+
+!
+!--          Write zu and zw data (vertical axes)
+             ALLOCATE( netcdf_data(0:nz+1) )
+
+             netcdf_data(0:nz+1) = zu(nzb:nzt+1)
+             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zu_xz(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ nz+2 /) )
+             CALL handle_netcdf_error( 'netcdf', 166 )
+
+             netcdf_data(0:nz+1) = zw(nzb:nzt+1)
+             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zw_xz(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ nz+2 /) )
+             CALL handle_netcdf_error( 'netcdf', 167 )
+
+             DEALLOCATE( netcdf_data )
+
+          ENDIF
 
 
@@ -2495,6 +2632,13 @@
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_time_xz(av), &
-                                            len = do2d_xz_time_count(av) )
+                                            len = ntime_count )
           CALL handle_netcdf_error( 'netcdf', 175 )
+
+!
+!--       For non-parallel output use the last output time level of the netcdf
+!--       file because the time dimension is unlimited. In case of parallel
+!--       output the variable ntime_count could get the value of 9*10E36 because
+!--       the time dimension is limited.
+          IF ( netcdf_data_format < 5 ) do2d_xz_time_count(av) = ntime_count
 
           nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_time_xz(av),    &
@@ -2517,4 +2661,20 @@
              RETURN
           ENDIF
+          
+          IF ( netcdf_data_format > 4 )  THEN
+             IF ( ntdim_2d_xz(av) > ntime_count )  THEN
+                message_string = 'netCDF file for cross sections ' // &
+                                 TRIM( var ) // ' from previous run found,' // &
+                                 '&but this file cannot be extended becaus' // &
+                                 'e the number of output time levels&has b' // &
+                                 'een increased compared to the previous s' // &
+                                 'imulation.' //                                &
+                                 '&New file is created instead.'
+                CALL message( 'define_netcdf_header', 'PA0390', 0, 1, 0, 6, 0 )
+                do2d_xz_time_count(av) = 0
+                extend = .FALSE.
+                RETURN
+             ENDIF
+          ENDIF
 
 !
@@ -2539,8 +2699,8 @@
                    ELSE
 !
-!--                   ATTENTION: Due to a probable bug in the netCDF4
-!--                              installation, independet output does not work
-!--                              A workaround is used in data_output_2d on those
-!--                              PEs having no data
+!--                   Test simulations showed that the output of cross sections
+!--                   by all PEs in data_output_2d using NF90_COLLECTIVE is 
+!--                   faster than the output by the first row of PEs in
+!--                   x-direction using NF90_INDEPENDENT.
                       nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
                                                      id_var_do2d(av,i), &
@@ -2610,8 +2770,16 @@
 
 !
-!--       Define time coordinate for yz sections (unlimited dimension)
-          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'time', NF90_UNLIMITED, &
-                                  id_dim_time_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 181 )
+!--       Define time coordinate for yz sections.
+!--       For parallel output the time dimensions has to be limited, otherwise
+!--       the performance drops significantly.
+          IF ( netcdf_data_format < 5 )  THEN
+             nc_stat = NF90_DEF_DIM( id_set_yz(av), 'time', NF90_UNLIMITED, &
+                                     id_dim_time_yz(av) )
+             CALL handle_netcdf_error( 'netcdf', 181 )
+          ELSE
+             nc_stat = NF90_DEF_DIM( id_set_yz(av), 'time', ntdim_2d_yz(av), &
+                                     id_dim_time_yz(av) )
+             CALL handle_netcdf_error( 'netcdf', 526 )
+          ENDIF
 
           nc_stat = NF90_DEF_VAR( id_set_yz(av), 'time', NF90_DOUBLE, &
@@ -2827,9 +2995,15 @@
                 CALL handle_netcdf_error( 'netcdf', 356 )
 #if defined( __netcdf4_parallel )
-!
-!--             Set independent io operations for parallel io. Collective io
-!--             is only allowed in case of a 1d-decomposition along y, because
-!--             otherwise, not all PEs have output data.
                 IF ( netcdf_data_format > 4 )  THEN
+!
+!--                Set no fill for every variable to increase performance.
+                   nc_stat = NF90_DEF_VAR_FILL( id_set_yz(av),     &
+                                                id_var_do2d(av,i), &
+                                                1, 0 )
+                   CALL handle_netcdf_error( 'netcdf', 535 )
+!
+!--                Set independent io operations for parallel io. Collective io
+!--                is only allowed in case of a 1d-decomposition along y,
+!--                because otherwise, not all PEs have output data.
                    IF ( npex == 1 )  THEN
                       nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av),     &
@@ -2838,8 +3012,8 @@
                    ELSE
 !
-!--                   ATTENTION: Due to a probable bug in the netCDF4
-!--                              installation, independet output does not work
-!--                              A workaround is used in data_output_2d on those
-!--                              PEs having no data
+!--                   Test simulations showed that the output of cross sections
+!--                   by all PEs in data_output_2d using NF90_COLLECTIVE is 
+!--                   faster than the output by the first row of PEs in
+!--                   y-direction using NF90_INDEPENDENT.
                       nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av),     &
                                                      id_var_do2d(av,i), &
@@ -2870,4 +3044,10 @@
 
 !
+!--       Set general no fill, otherwise the performance drops significantly for
+!--       parallel output.
+          nc_stat = NF90_SET_FILL( id_set_yz(av), NF90_NOFILL, oldmode )
+          CALL handle_netcdf_error( 'netcdf', 531 )
+
+!
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_yz(av) )
@@ -2875,83 +3055,93 @@
 
 !
-!--       Write axis data: x_yz, y, zu, zw
-          ALLOCATE( netcdf_data(1:ns) )
-
-!
-!--       Write x_yz data (shifted by +dx/2)
-          DO  i = 1, ns
-             IF( section(i,3) == -1 )  THEN
-                netcdf_data(i) = -1.0  ! section averaged along x
-             ELSE
-                netcdf_data(i) = ( section(i,3) + 0.5 ) * dx
-             ENDIF
-          ENDDO
-          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_x_yz(av), netcdf_data, &
-                                  start = (/ 1 /), count = (/ ns /) )
-          CALL handle_netcdf_error( 'netcdf', 202 ) 
-
-!
-!--       Write x_yz data (xu grid)
-          DO  i = 1, ns
-             IF( section(i,3) == -1 )  THEN
-                netcdf_data(i) = -1.0  ! section averaged along x
-             ELSE
-                netcdf_data(i) = section(i,3) * dx
-             ENDIF
-          ENDDO
-          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_xu_yz(av), netcdf_data, &
-                                  start = (/ 1 /), count = (/ ns /) )
-          CALL handle_netcdf_error( 'netcdf', 383 ) 
-
-!
-!--       Write gridpoint number data
-          netcdf_data(1:ns) = section(1:ns,3)
-          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_ind_x_yz(av), &
-                                  netcdf_data, start = (/ 1 /),       &
-                                  count = (/ ns /) )
-          CALL handle_netcdf_error( 'netcdf', 203 ) 
-
-          DEALLOCATE( netcdf_data )
-
-!
-!--       Write data for y (shifted by +dy/2) and yv axis
-          ALLOCATE( netcdf_data(0:ny+1) )
-
-          DO  j = 0, ny+1
-             netcdf_data(j) = ( j + 0.5 ) * dy
-          ENDDO
-
-          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_y_yz(av), netcdf_data, &
-                                  start = (/ 1 /), count = (/ ny+2 /) )
-           CALL handle_netcdf_error( 'netcdf', 204 ) 
-
-          DO  j = 0, ny+1
-             netcdf_data(j) = j * dy
-          ENDDO
-
-          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_yv_yz(av), &
-                                  netcdf_data, start = (/ 1 /),    &
-                                  count = (/ ny+2 /) )
-          CALL handle_netcdf_error( 'netcdf', 384 ) 
-
-          DEALLOCATE( netcdf_data )
-
-!
-!--       Write zu and zw data (vertical axes)
-          ALLOCATE( netcdf_data(0:nz+1) )
-
-          netcdf_data(0:nz+1) = zu(nzb:nzt+1)
-          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zu_yz(av), &
-                                  netcdf_data, start = (/ 1 /),    &
-                                  count = (/ nz+2 /) )
-          CALL handle_netcdf_error( 'netcdf', 205 ) 
-
-          netcdf_data(0:nz+1) = zw(nzb:nzt+1)
-          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zw_yz(av), &
-                                  netcdf_data, start = (/ 1 /),    &
-                                  count = (/ nz+2 /) )
-          CALL handle_netcdf_error( 'netcdf', 206 ) 
-
-          DEALLOCATE( netcdf_data )
+!--       These data are only written by PE0 for parallel output to increase
+!--       the performance.
+          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
+
+!
+!--          Write axis data: x_yz, y, zu, zw
+             ALLOCATE( netcdf_data(1:ns) )
+
+!
+!--          Write x_yz data (shifted by +dx/2)
+             DO  i = 1, ns
+                IF( section(i,3) == -1 )  THEN
+                   netcdf_data(i) = -1.0  ! section averaged along x
+                ELSE
+                   netcdf_data(i) = ( section(i,3) + 0.5 ) * dx
+                ENDIF
+             ENDDO
+             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_x_yz(av), &
+                                     netcdf_data, start = (/ 1 /),   &
+                                     count = (/ ns /) )
+             CALL handle_netcdf_error( 'netcdf', 202 ) 
+
+!
+!--          Write x_yz data (xu grid)
+             DO  i = 1, ns
+                IF( section(i,3) == -1 )  THEN
+                   netcdf_data(i) = -1.0  ! section averaged along x
+                ELSE
+                   netcdf_data(i) = section(i,3) * dx
+                ENDIF
+             ENDDO
+             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_xu_yz(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ ns /) )
+             CALL handle_netcdf_error( 'netcdf', 383 ) 
+
+!
+!--          Write gridpoint number data
+             netcdf_data(1:ns) = section(1:ns,3)
+             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_ind_x_yz(av), &
+                                     netcdf_data, start = (/ 1 /),       &
+                                     count = (/ ns /) )
+             CALL handle_netcdf_error( 'netcdf', 203 ) 
+
+             DEALLOCATE( netcdf_data )
+
+!
+!--          Write data for y (shifted by +dy/2) and yv axis
+             ALLOCATE( netcdf_data(0:ny+1) )
+
+             DO  j = 0, ny+1
+                netcdf_data(j) = ( j + 0.5 ) * dy
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_y_yz(av), &
+                                     netcdf_data, start = (/ 1 /),   &
+                                     count = (/ ny+2 /) )
+             CALL handle_netcdf_error( 'netcdf', 204 ) 
+
+             DO  j = 0, ny+1
+                netcdf_data(j) = j * dy
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_yv_yz(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ ny+2 /) )
+             CALL handle_netcdf_error( 'netcdf', 384 ) 
+
+             DEALLOCATE( netcdf_data )
+
+!
+!--          Write zu and zw data (vertical axes)
+             ALLOCATE( netcdf_data(0:nz+1) )
+
+             netcdf_data(0:nz+1) = zu(nzb:nzt+1)
+             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zu_yz(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ nz+2 /) )
+             CALL handle_netcdf_error( 'netcdf', 205 ) 
+
+             netcdf_data(0:nz+1) = zw(nzb:nzt+1)
+             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zw_yz(av), &
+                                     netcdf_data, start = (/ 1 /),    &
+                                     count = (/ nz+2 /) )
+             CALL handle_netcdf_error( 'netcdf', 206 ) 
+
+             DEALLOCATE( netcdf_data )
+
+          ENDIF
 
 
@@ -3080,6 +3270,13 @@
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_yz(av), id_dim_time_yz(av), &
-                                            len = do2d_yz_time_count(av) )
+                                            len = ntime_count )
           CALL handle_netcdf_error( 'netcdf', 214 )
+
+!
+!--       For non-parallel output use the last output time level of the netcdf
+!--       file because the time dimension is unlimited. In case of parallel
+!--       output the variable ntime_count could get the value of 9*10E36 because
+!--       the time dimension is limited.
+          IF ( netcdf_data_format < 5 ) do2d_yz_time_count(av) = ntime_count
 
           nc_stat = NF90_GET_VAR( id_set_yz(av), id_var_time_yz(av),    &
@@ -3103,4 +3300,20 @@
           ENDIF
 
+          IF ( netcdf_data_format > 4 )  THEN
+             IF ( ntdim_2d_yz(av) > ntime_count )  THEN
+                message_string = 'netCDF file for cross sections ' //          &
+                                 TRIM( var ) // ' from previous run found,' // &
+                                 '&but this file cannot be extended becaus' // &
+                                 'e the number of output time levels&has b' // &
+                                 'een increased compared to the previous s' // &
+                                 'imulation.' //                               &
+                                 '&New file is created instead.'
+                CALL message( 'define_netcdf_header', 'PA0391', 0, 1, 0, 6, 0 )
+                do2d_yz_time_count(av) = 0
+                extend = .FALSE.
+                RETURN
+             ENDIF
+          ENDIF
+
 !
 !--       Dataset seems to be extendable.
@@ -3124,8 +3337,8 @@
                    ELSE
 !
-!--                   ATTENTION: Due to a probable bug in the netCDF4
-!--                              installation, independet output does not work
-!--                              A workaround is used in data_output_2d on those
-!--                              PEs having no data
+!--                   Test simulations showed that the output of cross sections
+!--                   by all PEs in data_output_2d using NF90_COLLECTIVE is 
+!--                   faster than the output by the first row of PEs in
+!--                   y-direction using NF90_INDEPENDENT.
                       nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av),     &
                                                      id_var_do2d(av,i), &
@@ -4651,5 +4864,5 @@
 !--    64bit-offset format
        nc_stat = NF90_CREATE( filename,                                        &
-                              OR( NF90_NOCLOBBER, NF90_64BIT_OFFSET ), id )
+                              IOR( NF90_NOCLOBBER, NF90_64BIT_OFFSET ), id )
 
 #if defined( __netcdf4 )
@@ -4658,5 +4871,5 @@
 !
 !--    netCDF4/HDF5 format
-       nc_stat = NF90_CREATE( filename, OR( NF90_NOCLOBBER, NF90_NETCDF4 ), id )
+       nc_stat = NF90_CREATE( filename, IOR( NF90_NOCLOBBER, NF90_NETCDF4 ), id )
 
     ELSEIF ( netcdf_data_format == 4  .OR.                                     &
@@ -4665,6 +4878,6 @@
 !--    netCDF4/HDF5 format with classic model flag
        nc_stat = NF90_CREATE( filename,                                        &
-                              OR( NF90_NOCLOBBER,                              &
-                              OR( NF90_CLASSIC_MODEL, NF90_HDF5 ) ), id )
+                              IOR( NF90_NOCLOBBER,                              &
+                              IOR( NF90_CLASSIC_MODEL, NF90_HDF5 ) ), id )
 
 #if defined( __netcdf4_parallel )
@@ -4672,5 +4885,7 @@
 !
 !--    netCDF4/HDF5 format, parallel
-       nc_stat = NF90_CREATE( filename, OR( NF90_NOCLOBBER, NF90_NETCDF4 ),    &
+       nc_stat = NF90_CREATE( filename,                                        &
+                              IOR( NF90_NOCLOBBER,                              &
+                              IOR( NF90_NETCDF4, NF90_MPIIO ) ),                &
                               id, COMM = comm2d, INFO = MPI_INFO_NULL )
 
@@ -4679,6 +4894,7 @@
 !--    netCDF4/HDF5 format with classic model flag, parallel
        nc_stat = NF90_CREATE( filename,                                        &
-                              OR( NF90_NOCLOBBER,                              &
-                              OR( NF90_CLASSIC_MODEL, NF90_HDF5 ) ),           &
+                              IOR( NF90_NOCLOBBER,                              &
+                              IOR( NF90_MPIIO,                                  &
+                              IOR( NF90_CLASSIC_MODEL, NF90_HDF5 ) ) ),         &
                               id, COMM = comm2d, INFO = MPI_INFO_NULL )
 
@@ -4720,6 +4936,6 @@
 #if defined( __netcdf4_parallel )
     ELSEIF ( netcdf_data_format > 4  .AND.  parallel )  THEN
-       nc_stat = NF90_OPEN( filename, NF90_WRITE, id, COMM = comm2d, &
-                            INFO = MPI_INFO_NULL )
+       nc_stat = NF90_OPEN( filename, IOR( NF90_WRITE, NF90_MPIIO ), id,  &
+                            COMM = comm2d, INFO = MPI_INFO_NULL )
 #endif
 #endif
Index: palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- palm/trunk/SOURCE/read_var_list.f90	(revision 1307)
+++ palm/trunk/SOURCE/read_var_list.f90	(revision 1308)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-!
+! +do2d_xy_time_count, do2d_xz_time_count, do2d_yz_time_count,
+! +do3d_time_count
 !
 ! Former revisions:
@@ -202,5 +203,5 @@
 !-- Make version number check first
     READ ( 13 )  version_on_file
-    binary_version = '3.9a'
+    binary_version = '3.9b'
     IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
        WRITE( message_string, * ) 'version mismatch concerning control ', &
@@ -383,4 +384,12 @@
           CASE ( 'dissipation_1d' )
              READ ( 13 )  dissipation_1d
+          CASE ( 'do2d_xy_time_count' )
+             READ ( 13 )  do2d_xy_time_count
+          CASE ( 'do2d_xz_time_count' )
+             READ ( 13 )  do2d_xz_time_count
+          CASE ( 'do2d_yz_time_count' )
+             READ ( 13 )  do2d_yz_time_count
+          CASE ( 'do3d_time_count' )
+             READ ( 13 )  do3d_time_count
           CASE ( 'dp_external' )
              READ ( 13 )  dp_external
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 1307)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 1308)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! ------------------
-! 
+! +netcdf_data_format_save
+! For masked data, parallel netcdf output is not tested so far, hence 
+! netcdf_data_format is switched back to non-paralell output.
 !
 ! Former revisions:
@@ -186,4 +188,5 @@
 
     CHARACTER (LEN=9) ::  time_to_string
+    INTEGER           ::  netcdf_data_format_save
 
 !
@@ -810,4 +813,9 @@
 !
 !--    masked data output
+!--    Parallel netcdf output is not tested so far for masked data, hence 
+!--    netcdf_data_format is switched back to non-paralell output.
+       netcdf_data_format_save = netcdf_data_format
+       IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
+       IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
        DO  mid = 1, masks
           IF ( time_domask(mid) >= dt_domask(mid) )  THEN
@@ -817,4 +825,5 @@
           ENDIF
        ENDDO
+       netcdf_data_format = netcdf_data_format_save
 
 !
@@ -826,7 +835,13 @@
           CALL data_output_2d( 'yz', 1 )
           CALL data_output_3d( 1 )
+!--       Parallel netcdf output is not tested so far for masked data, hence 
+!--       netcdf_data_format is switched back to non-paralell output.
+          netcdf_data_format_save = netcdf_data_format
+          IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
+          IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
           DO  mid = 1, masks
              CALL data_output_mask( 1 )
           ENDDO
+          netcdf_data_format = netcdf_data_format_save
           time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
        ENDIF
Index: palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- palm/trunk/SOURCE/write_var_list.f90	(revision 1307)
+++ palm/trunk/SOURCE/write_var_list.f90	(revision 1308)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-!
+! +do2d_xy_time_count, do2d_xz_time_count, do2d_yz_time_count,
+! +do3d_time_count
 !
 ! Former revisions:
@@ -176,5 +177,5 @@
 
 
-    binary_version = '3.9a'
+    binary_version = '3.9b'
 
     WRITE ( 14 )  binary_version
@@ -297,4 +298,12 @@
     WRITE ( 14 )  'dissipation_1d                '
     WRITE ( 14 )  dissipation_1d
+    WRITE ( 14 )  'do2d_xy_time_count            '
+    WRITE ( 14 )  do2d_xy_time_count
+    WRITE ( 14 )  'do2d_xz_time_count            '
+    WRITE ( 14 )  do2d_xz_time_count
+    WRITE ( 14 )  'do2d_yz_time_count            '
+    WRITE ( 14 )  do2d_yz_time_count
+    WRITE ( 14 )  'do3d_time_count               '
+    WRITE ( 14 )  do3d_time_count
     WRITE ( 14 )  'dp_external                   '
     WRITE ( 14 )  dp_external
