Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 4069)
+++ palm/trunk/SOURCE/Makefile	(revision 4070)
@@ -25,6 +25,9 @@
 # -----------------
 # $Id$
+# Add new data output modules
+#
+# 4048 2019-06-21 21:00:21Z knoop
 # Moved turbulence_closure_mod calls into module_interface
-# 
+#
 # 4047 2019-06-21 18:58:09Z knoop
 # Introduction of the dynamics module and dissolution of swap_timelevel
@@ -597,4 +600,8 @@
 	cuda_fft_interfaces.f90 \
 	data_log.f90 \
+	data_output_binary_module.f90 \
+	data_output_module.f90 \
+	data_output_netcdf4_parallel_module.f90 \
+	data_output_netcdf4_serial_module.f90 \
 	data_output_2d.f90 \
 	data_output_3d.f90 \
@@ -903,4 +910,11 @@
 	mod_kinds.o \
 	modules.o
+data_output_binary_module.o: \
+	mod_kinds.o
+data_output_module.o: \
+	data_output_binary_module.o \
+	data_output_netcdf4_parallel_module.o \
+	data_output_netcdf4_serial_module.o \
+	mod_kinds.o
 data_output_mask.o: \
 	basic_constants_and_equations_mod.o \
@@ -915,4 +929,8 @@
 	salsa_mod.o \
 	surface_mod.o
+data_output_netcdf4_parallel_module.o: \
+	mod_kinds.o
+data_output_netcdf4_serial_module.o: \
+	mod_kinds.o
 data_output_profiles.o: \
 	cpulog_mod.o \
Index: palm/trunk/SOURCE/data_output_binary_module.f90
===================================================================
--- palm/trunk/SOURCE/data_output_binary_module.f90	(revision 4070)
+++ palm/trunk/SOURCE/data_output_binary_module.f90	(revision 4070)
@@ -0,0 +1,734 @@
+!> @file data_output_binary_module.f90
+!--------------------------------------------------------------------------------------------------!
+! This file is part of the PALM model system.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 2019-2019 Leibniz Universitaet Hannover
+!--------------------------------------------------------------------------------------------------!
+!
+! Current revisions:
+! ------------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
+! Initial revision
+!
+!
+! Authors:
+! --------
+!> @author: Tobias Gronemeier
+!
+! Description:
+! ------------
+!> Binary output module to write output data into binary files.
+!>
+!> @todo Think of removing 'is_init' as its value can be derived from the return
+!>       value of 'return_value'.
+!> @todo Get return value of write statements.
+!--------------------------------------------------------------------------------------------------!
+MODULE data_output_binary_module
+
+   USE kinds
+
+#if defined( __parallel )
+#if defined( __mpifh )
+   INCLUDE "mpif.h"
+#else
+   USE MPI
+#endif
+#endif
+
+   IMPLICIT NONE
+
+   INTEGER(iwp), PARAMETER ::  charlen = 100_iwp  !< maximum length of character variables
+
+   CHARACTER(LEN=*), PARAMETER ::  config_file_name = 'BINARY_TO_NETCDF_CONFIG'  !< name of config file
+   CHARACTER(LEN=*), PARAMETER ::  prefix = 'BIN_'                               !< file prefix for binary files
+
+   CHARACTER(LEN=800)          ::  internal_error_message = ''  !< string containing the last error message
+   CHARACTER(LEN=800)          ::  temp_string                  !< dummy string
+
+   INTEGER(iwp) ::  binary_file_lowest_unit = 1000  !< lowest unit number of all binary files created by this module
+   INTEGER(iwp) ::  config_file_unit                !< unit number of config file
+   INTEGER(iwp) ::  debug_output_unit               !< Fortran Unit Number of the debug-output file
+   INTEGER(iwp) ::  global_id_in_file = -1          !< value of global ID within a file
+   INTEGER(iwp) ::  next_available_unit             !< next unit number available for new file
+
+   INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  files_highest_var_id  !< highest assigned ID of variable or dimension in a file
+
+   LOGICAL ::  binary_open_file_first_call = .TRUE.  !< true if binary_open_file routine was not called yet
+   LOGICAL ::  config_file_open = .FALSE.            !< true if config file is opened and not closed
+   LOGICAL ::  print_debug_output = .FALSE.          !< if true, debug output is printed
+
+   SAVE
+
+   PRIVATE
+
+   INTERFACE binary_init_module
+      MODULE PROCEDURE binary_init_module
+   END INTERFACE binary_init_module
+
+   INTERFACE binary_open_file
+      MODULE PROCEDURE binary_open_file
+   END INTERFACE binary_open_file
+
+   INTERFACE binary_init_dimension
+      MODULE PROCEDURE binary_init_dimension
+   END INTERFACE binary_init_dimension
+
+   INTERFACE binary_init_variable
+      MODULE PROCEDURE binary_init_variable
+   END INTERFACE binary_init_variable
+
+   INTERFACE binary_write_attribute
+      MODULE PROCEDURE binary_write_attribute
+   END INTERFACE binary_write_attribute
+
+   INTERFACE binary_init_end
+      MODULE PROCEDURE binary_init_end
+   END INTERFACE binary_init_end
+
+   INTERFACE binary_write_variable
+      MODULE PROCEDURE binary_write_variable
+   END INTERFACE binary_write_variable
+
+   INTERFACE binary_finalize
+      MODULE PROCEDURE binary_finalize
+   END INTERFACE binary_finalize
+
+   INTERFACE binary_get_error_message
+      MODULE PROCEDURE binary_get_error_message
+   END INTERFACE binary_get_error_message
+
+   PUBLIC &
+      binary_finalize, &
+      binary_get_error_message, &
+      binary_init_dimension, &
+      binary_init_end, &
+      binary_init_module, &
+      binary_init_variable, &
+      binary_open_file, &
+      binary_write_attribute, &
+      binary_write_variable
+
+
+CONTAINS
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize data-output module.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE binary_init_module( program_debug_output_unit, debug_output, dom_global_id )
+
+   INTEGER(iwp), INTENT(IN) ::  dom_global_id              !< global id within a file defined by DOM
+   INTEGER(iwp), INTENT(IN) ::  program_debug_output_unit  !< file unit number for debug output
+
+   LOGICAL, INTENT(IN) ::  debug_output  !< if true, debug output is printed
+
+
+   debug_output_unit = program_debug_output_unit
+
+   print_debug_output = debug_output
+
+   global_id_in_file = dom_global_id
+
+
+END SUBROUTINE binary_init_module
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Open binary file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE binary_open_file( filename, file_id, return_value )
+
+   CHARACTER(LEN=charlen)             ::  bin_filename = ''  !< actual name of binary file
+   CHARACTER(LEN=charlen), INTENT(IN) ::  filename           !< name of file
+   CHARACTER(LEN=7)                   ::  myid_char          !< string containing value of myid with leading zeros
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'binary_open_file'  !< name of this routine
+
+   INTEGER(iwp), INTENT(OUT) ::  file_id       !< file ID
+   INTEGER(iwp)              ::  myid          !< id of local processor id (MPI rank)
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+   INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  files_highest_var_id_tmp  !< temporary list of given variable IDs in file
+
+   LOGICAL ::  file_exists  !< true if file to be opened already exists
+
+
+   return_value = 0
+
+#if defined( __parallel )
+   CALL MPI_COMM_RANK( MPI_COMM_WORLD, myid, return_value )
+   IF ( return_value == 0 )  THEN
+      WRITE( myid_char, '("_",I6.6)' )  myid
+   ELSE
+      CALL internal_message( 'debug', routine_name // 'MPI_COMM_RANK error' )
+   ENDIF
+#else
+   myid = 0
+   myid_char = '_' // REPEAT('0', 6)
+#endif
+
+   !-- Open binary config file for combining script
+   IF ( return_value == 0  .AND.  binary_open_file_first_call )  THEN
+
+      binary_open_file_first_call = .FALSE.
+      config_file_unit = binary_file_lowest_unit
+
+      IF ( myid == 0 )  THEN
+
+         !-- Remove any pre-existing file
+         INQUIRE( FILE=TRIM( config_file_name ), EXIST=file_exists )
+
+         IF ( file_exists )  THEN
+            CALL internal_message( 'debug', routine_name //              &
+                                            ': Remove existing file ' // &
+                                            TRIM( config_file_name ) )
+            CALL SYSTEM( 'rm ' // TRIM( config_file_name ) )
+         ENDIF
+
+         OPEN( config_file_unit, FILE=TRIM( config_file_name ), &
+               FORM='UNFORMATTED', STATUS='NEW', IOSTAT=return_value )
+
+         IF ( return_value == 0 )  THEN
+
+            config_file_open = .TRUE.
+
+            !-- Write some general information to config file
+            WRITE( config_file_unit )  LEN( prefix )
+            WRITE( config_file_unit )  prefix
+            WRITE( config_file_unit )  charlen
+            WRITE( config_file_unit )  global_id_in_file
+
+         ELSE
+
+            return_value = 1
+            CALL internal_message( 'error', routine_name // ': could not create config' )
+
+         ENDIF
+
+      ENDIF
+
+      next_available_unit = binary_file_lowest_unit + 1
+
+   ENDIF
+
+   !-- Initialize output file: open, write header, initialize variable/dimension IDs
+   IF ( return_value == 0 )  THEN
+
+      bin_filename = prefix // TRIM( filename ) // myid_char
+
+      !-- Remove any pre-existing file
+      INQUIRE( FILE=TRIM( bin_filename ), EXIST=file_exists )
+
+      IF ( file_exists )  THEN
+         CALL internal_message( 'debug', routine_name // &
+                                         ': remove existing file ' // TRIM( bin_filename ) )
+         CALL SYSTEM( 'rm ' // TRIM( bin_filename ) )
+      ENDIF
+
+      !-- Open binary file
+      CALL internal_message( 'debug', routine_name // ': open file ' // TRIM( bin_filename ) )
+      OPEN ( next_available_unit, FILE=TRIM( bin_filename ), &
+             FORM='UNFORMATTED', STATUS='NEW', IOSTAT=return_value )
+
+      IF ( return_value == 0 )  THEN
+
+         !-- Add filename to config file
+         IF ( myid == 0 )  THEN
+            WRITE( config_file_unit )  filename
+         ENDIF
+
+         !-- Save file ID and increase next file unit number
+         file_id = next_available_unit
+         next_available_unit = next_available_unit + 1
+
+         !-- Write some meta data to file
+         WRITE ( file_id )  charlen
+         WRITE ( file_id )  file_id
+         WRITE ( file_id )  filename
+
+         !-- Extend file-variable/dimension-ID list by 1 and set it to 0 for new file.
+         IF ( ALLOCATED( files_highest_var_id ) )  THEN
+            ALLOCATE( files_highest_var_id_tmp(SIZE( files_highest_var_id )) )
+            files_highest_var_id_tmp = files_highest_var_id
+            DEALLOCATE( files_highest_var_id )
+            ALLOCATE( files_highest_var_id(binary_file_lowest_unit+1:file_id) )
+            files_highest_var_id(:file_id-1) = files_highest_var_id_tmp
+            DEALLOCATE( files_highest_var_id_tmp )
+         ELSE
+            ALLOCATE( files_highest_var_id(binary_file_lowest_unit+1:file_id) )
+         ENDIF
+         files_highest_var_id(file_id) = 0_iwp
+
+      ELSE
+         return_value = 1
+         CALL internal_message( 'error', routine_name // &
+                                         ': could not open file "' // TRIM( filename ) // '"')
+      ENDIF
+
+   ENDIF
+
+END SUBROUTINE binary_open_file
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Write attribute to file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE binary_write_attribute( file_id, var_id, att_name, att_value_char, &
+              att_value_int8, att_value_int16, att_value_int32,               &
+              att_value_real32, att_value_real64, return_value )
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'binary_write_attribute'  !< name of this routine
+
+   CHARACTER(LEN=charlen), INTENT(IN)           ::  att_name        !< name of attribute
+   CHARACTER(LEN=charlen), INTENT(IN), OPTIONAL ::  att_value_char  !< value of attribute
+   CHARACTER(LEN=charlen)                       ::  att_type        !< data type of attribute
+   CHARACTER(LEN=charlen)                       ::  out_str         !< output string
+
+   INTEGER(KIND=1), INTENT(IN), OPTIONAL ::  att_value_int8   !< value of attribute
+   INTEGER(KIND=2), INTENT(IN), OPTIONAL ::  att_value_int16  !< value of attribute
+   INTEGER(KIND=4), INTENT(IN), OPTIONAL ::  att_value_int32  !< value of attribute
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp), INTENT(IN)  ::  var_id        !< variable ID
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+   REAL(KIND=4), INTENT(IN), OPTIONAL ::  att_value_real32  !< value of attribute
+   REAL(KIND=8), INTENT(IN), OPTIONAL ::  att_value_real64  !< value of attribute
+
+
+   return_value = 0
+
+   CALL internal_message( 'debug', TRIM( routine_name ) // &
+                                   ': write attribute ' // TRIM( att_name ) )
+
+   !-- Write attribute to file
+   out_str = 'attribute'
+   WRITE( file_id )  out_str
+
+   WRITE( file_id )  var_id
+   WRITE( file_id )  att_name
+
+   IF ( PRESENT( att_value_char ) )  THEN
+      att_type = 'char'
+      WRITE( file_id )  att_type
+      WRITE( file_id )  att_value_char
+   ELSEIF ( PRESENT( att_value_int8 ) )  THEN
+      att_type = 'int8'
+      WRITE( file_id )  att_type
+      WRITE( file_id )  att_value_int8
+   ELSEIF ( PRESENT( att_value_int16 ) )  THEN
+      att_type = 'int16'
+      WRITE( file_id )  att_type
+      WRITE( file_id )  att_value_int16
+   ELSEIF ( PRESENT( att_value_int32 ) )  THEN
+      att_type = 'int32'
+      WRITE( file_id )  att_type
+      WRITE( file_id )  att_value_int32
+   ELSEIF ( PRESENT( att_value_real32 ) )  THEN
+      att_type = 'real32'
+      WRITE( file_id )  att_type
+      WRITE( file_id )  att_value_real32
+   ELSEIF ( PRESENT( att_value_real64 ) )  THEN
+      att_type = 'real64'
+      WRITE( file_id )  att_type
+      WRITE( file_id )  att_value_real64
+   ELSE
+      return_value = 1
+      CALL internal_message( 'error', TRIM( routine_name ) // &
+                                      ': attribute "' // TRIM( att_name ) // '": no value given' )
+   ENDIF
+
+END SUBROUTINE binary_write_attribute
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize dimension. Write information in file header and save dimension
+!> values to be later written to file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE binary_init_dimension( file_id, dim_id, var_id, &
+              dim_name, dim_type, dim_length, is_init, return_value )
+
+   CHARACTER(LEN=charlen), INTENT(IN) ::  dim_name  !< name of dimension
+   CHARACTER(LEN=charlen), INTENT(IN) ::  dim_type  !< data type of dimension
+   CHARACTER(LEN=charlen)             ::  out_str   !< output string
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'binary_init_dimension'  !< name of this routine
+
+   INTEGER(iwp), INTENT(OUT) ::  dim_id        !< dimension ID
+   INTEGER(iwp), INTENT(IN)  ::  dim_length    !< length of dimension
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+   INTEGER(iwp), INTENT(OUT) ::  var_id        !< variable ID
+
+   LOGICAL, INTENT(OUT) ::  is_init  !< true if dimension is initialized
+
+
+   return_value = 0
+
+   CALL internal_message( 'debug', routine_name // ': init dimension ' // TRIM( dim_name ) )
+
+   !-- Assign dimension ID
+   dim_id = files_highest_var_id( file_id ) + 1
+   files_highest_var_id( file_id ) = dim_id
+
+   !-- Define dimension in file
+   out_str = 'dimension'
+   WRITE( file_id )  out_str
+   WRITE( file_id )  dim_name
+   WRITE( file_id )  dim_id
+   WRITE( file_id )  dim_type
+   WRITE( file_id )  dim_length
+
+   !-- Define variable associated with dimension
+   CALL binary_init_variable( file_id, var_id, dim_name, dim_type, &
+                              (/dim_id/), is_init, .TRUE., return_value )
+   IF ( return_value /= 0 )  THEN
+      is_init = .FALSE.
+      CALL internal_message( 'error', routine_name // &
+                                      ': init dimension "' // TRIM( dim_name ) // '"' )
+   ENDIF
+
+END SUBROUTINE binary_init_dimension
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize variable. Write information of variable into file header.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE binary_init_variable( file_id, var_id, var_name, var_type, &
+                                 var_dim_ids, is_init, is_global, return_value )
+
+   CHARACTER(LEN=charlen)             ::  out_str   !< output string
+   CHARACTER(LEN=charlen), INTENT(IN) ::  var_name  !< name of variable
+   CHARACTER(LEN=charlen), INTENT(IN) ::  var_type  !< data type of variable
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'binary_init_variable'  !< name of this routine
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp), INTENT(OUT) ::  var_id        !< variable ID
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+   INTEGER(iwp), DIMENSION(:), INTENT(IN) ::  var_dim_ids  !< list of dimension IDs used by variable
+
+   LOGICAL, INTENT(IN)  ::  is_global  !< true if variable is global (same on all PE)
+   LOGICAL, INTENT(OUT) ::  is_init    !< true if variable is initialized
+
+
+   return_value = 0
+
+   CALL internal_message( 'debug', routine_name // ': init variable ' // TRIM( var_name ) )
+
+   IF ( is_global )  CONTINUE  ! reqired to prevent compiler warning
+
+   !-- Assign variable ID
+   var_id = files_highest_var_id( file_id ) + 1
+   files_highest_var_id( file_id ) = var_id
+
+   !-- Write variable information in file
+   out_str = 'variable'
+   WRITE( file_id )  out_str
+   WRITE( file_id )  var_name
+   WRITE( file_id )  var_id
+   WRITE( file_id )  var_type
+   WRITE( file_id )  SIZE( var_dim_ids )
+   WRITE( file_id )  var_dim_ids
+
+   !-- Variable is initialised
+   is_init = return_value == 0
+
+END SUBROUTINE binary_init_variable
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Leave file definition state.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE binary_init_end( file_id, return_value )
+
+   CHARACTER(LEN=charlen) ::  out_str  !< output string
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'binary_init_end'  !< name of this routine
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+
+   return_value = 0
+
+   WRITE( temp_string, * ) file_id
+   CALL internal_message( 'debug', &
+                          routine_name // &
+                          ': finalize file definition (file_id=' // TRIM( temp_string ) // ')' )
+
+   out_str = '*** end file header ***'
+   WRITE( file_id )  out_str
+
+END SUBROUTINE binary_init_end
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Write variable to file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE binary_write_variable(                                         &
+              file_id, var_id, bounds_start, bounds_end, bounds_origin,   &
+              do_output, is_global,                                       &
+              var_int8_0d,   var_int8_1d,   var_int8_2d,   var_int8_3d,   &
+              var_int16_0d,  var_int16_1d,  var_int16_2d,  var_int16_3d,  &
+              var_int32_0d,  var_int32_1d,  var_int32_2d,  var_int32_3d,  &
+              var_intwp_0d,  var_intwp_1d,  var_intwp_2d,  var_intwp_3d,  &
+              var_real32_0d, var_real32_1d, var_real32_2d, var_real32_3d, &
+              var_real64_0d, var_real64_1d, var_real64_2d, var_real64_3d, &
+              var_realwp_0d, var_realwp_1d, var_realwp_2d, var_realwp_3d, &
+              return_value )
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'binary_write_variable'  !< name of this routine
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+   INTEGER(iwp), INTENT(IN)  ::  var_id        !< variable ID
+
+   INTEGER(iwp), DIMENSION(:), INTENT(IN) ::  bounds_origin  !< starting index of each dimension
+   INTEGER(iwp), DIMENSION(:), INTENT(IN) ::  bounds_end     !< ending index of variable
+   INTEGER(iwp), DIMENSION(:), INTENT(IN) ::  bounds_start   !< starting index of variable
+
+   INTEGER(KIND=1), POINTER, INTENT(IN), OPTIONAL                               ::  var_int8_0d  !< output variable
+   INTEGER(KIND=1), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_int8_1d  !< output variable
+   INTEGER(KIND=1), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_int8_2d  !< output variable
+   INTEGER(KIND=1), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_int8_3d  !< output variable
+
+   INTEGER(KIND=2), POINTER, INTENT(IN), OPTIONAL                               ::  var_int16_0d  !< output variable
+   INTEGER(KIND=2), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_int16_1d  !< output variable
+   INTEGER(KIND=2), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_int16_2d  !< output variable
+   INTEGER(KIND=2), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_int16_3d  !< output variable
+
+   INTEGER(KIND=4), POINTER, INTENT(IN), OPTIONAL                               ::  var_int32_0d  !< output variable
+   INTEGER(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_int32_1d  !< output variable
+   INTEGER(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_int32_2d  !< output variable
+   INTEGER(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_int32_3d  !< output variable
+
+   INTEGER(iwp), POINTER, INTENT(IN), OPTIONAL                               ::  var_intwp_0d  !< output variable
+   INTEGER(iwp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_intwp_1d  !< output variable
+   INTEGER(iwp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_intwp_2d  !< output variable
+   INTEGER(iwp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_intwp_3d  !< output variable
+
+   LOGICAL, INTENT(IN) ::  do_output  !< write output only if do_output = true
+   LOGICAL, INTENT(IN) ::  is_global  !< true if variable is global (same on all PE)
+
+   REAL(KIND=4), POINTER, INTENT(IN), OPTIONAL                               ::  var_real32_0d  !< output variable
+   REAL(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_real32_1d  !< output variable
+   REAL(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_real32_2d  !< output variable
+   REAL(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_real32_3d  !< output variable
+
+   REAL(KIND=8), POINTER, INTENT(IN), OPTIONAL                               ::  var_real64_0d  !< output variable
+   REAL(KIND=8), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_real64_1d  !< output variable
+   REAL(KIND=8), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_real64_2d  !< output variable
+   REAL(KIND=8), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_real64_3d  !< output variable
+
+   REAL(wp), POINTER, INTENT(IN), OPTIONAL                               ::  var_realwp_0d  !< output variable
+   REAL(wp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_realwp_1d  !< output variable
+   REAL(wp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_realwp_2d  !< output variable
+   REAL(wp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_realwp_3d  !< output variable
+
+
+   return_value = 0
+
+   WRITE( temp_string, '(": write variable ",I6," into file ",I6)' ) var_id, file_id
+   CALL internal_message( 'debug', routine_name // TRIM( temp_string ) )
+
+   IF ( is_global )  CONTINUE  ! reqired to prevent compiler warning
+
+   IF ( do_output )  THEN
+      WRITE( file_id )  var_id
+      WRITE( file_id )  bounds_start
+      WRITE( file_id )  bounds_end
+      WRITE( file_id )  bounds_origin
+      !-- 8bit integer output
+      IF ( PRESENT( var_int8_0d ) )  THEN
+         WRITE( file_id )  var_int8_0d
+      ELSEIF ( PRESENT( var_int8_1d ) )  THEN
+         WRITE( file_id )  var_int8_1d
+      ELSEIF ( PRESENT( var_int8_2d ) )  THEN
+         WRITE( file_id )  var_int8_2d
+      ELSEIF ( PRESENT( var_int8_3d ) )  THEN
+         WRITE( file_id )  var_int8_3d
+      !-- 16bit integer output
+      ELSEIF ( PRESENT( var_int16_0d ) )  THEN
+         WRITE( file_id )  var_int16_0d
+      ELSEIF ( PRESENT( var_int16_1d ) )  THEN
+         WRITE( file_id )  var_int16_1d
+      ELSEIF ( PRESENT( var_int16_2d ) )  THEN
+         WRITE( file_id )  var_int16_2d
+      ELSEIF ( PRESENT( var_int16_3d ) )  THEN
+         WRITE( file_id )  var_int16_3d
+      !-- 32bit integer output
+      ELSEIF ( PRESENT( var_int32_0d ) )  THEN
+         WRITE( file_id )  var_int32_0d
+      ELSEIF ( PRESENT( var_int32_1d ) )  THEN
+         WRITE( file_id )  var_int32_1d
+      ELSEIF ( PRESENT( var_int32_2d ) )  THEN
+         WRITE( file_id )  var_int32_2d
+      ELSEIF ( PRESENT( var_int32_3d ) )  THEN
+         WRITE( file_id )  var_int32_3d
+      !-- working-precision integer output
+      ELSEIF ( PRESENT( var_intwp_0d ) )  THEN
+         WRITE( file_id )  var_intwp_0d
+      ELSEIF ( PRESENT( var_intwp_1d ) )  THEN
+         WRITE( file_id )  var_intwp_1d
+      ELSEIF ( PRESENT( var_intwp_2d ) )  THEN
+         WRITE( file_id )  var_intwp_2d
+      ELSEIF ( PRESENT( var_intwp_3d ) )  THEN
+         WRITE( file_id )  var_intwp_3d
+      !-- 32bit real output
+      ELSEIF ( PRESENT( var_real32_0d ) )  THEN
+         WRITE( file_id )  var_real32_0d
+      ELSEIF ( PRESENT( var_real32_1d ) )  THEN
+         WRITE( file_id )  var_real32_1d
+      ELSEIF ( PRESENT( var_real32_2d ) )  THEN
+         WRITE( file_id )  var_real32_2d
+      ELSEIF ( PRESENT( var_real32_3d ) )  THEN
+         WRITE( file_id )  var_real32_3d
+      !-- 64bit real output
+      ELSEIF ( PRESENT( var_real64_0d ) )  THEN
+         WRITE( file_id )  var_real64_0d
+      ELSEIF ( PRESENT( var_real64_1d ) )  THEN
+         WRITE( file_id )  var_real64_1d
+      ELSEIF ( PRESENT( var_real64_2d ) )  THEN
+         WRITE( file_id )  var_real64_2d
+      ELSEIF ( PRESENT( var_real64_3d ) )  THEN
+         WRITE( file_id )  var_real64_3d
+      !-- working-precision real output
+      ELSEIF ( PRESENT( var_realwp_0d ) )  THEN
+         WRITE( file_id )  var_realwp_0d
+      ELSEIF ( PRESENT( var_realwp_1d ) )  THEN
+         WRITE( file_id )  var_realwp_1d
+      ELSEIF ( PRESENT( var_realwp_2d ) )  THEN
+         WRITE( file_id )  var_realwp_2d
+      ELSEIF ( PRESENT( var_realwp_3d ) )  THEN
+         WRITE( file_id )  var_realwp_3d
+      ELSE
+         return_value = 1
+         WRITE( temp_string, '(": variable_id=",I6,", file_id=",I6)' ) var_id, file_id
+         CALL internal_message( 'error', routine_name // &
+                                         TRIM( temp_string ) // ': no values given' )
+      ENDIF
+
+   ENDIF
+
+END SUBROUTINE binary_write_variable
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Close opened files.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE binary_finalize( file_id, return_value )
+
+   CHARACTER(LEN=charlen) ::  out_str  !< output string
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'binary_finalize'  !< name of this routine
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+
+   IF ( config_file_open )  THEN
+
+      out_str = '*** end config file ***'
+      WRITE( config_file_unit )  out_str
+
+      CLOSE( config_file_unit, IOSTAT=return_value )
+
+      IF ( return_value /= 0 )  THEN
+         CALL internal_message( 'error', routine_name // ': cannot close configuration file' )
+      ELSE
+         config_file_open = .FALSE.
+      ENDIF
+
+   ELSE
+
+      return_value = 0
+
+   ENDIF
+
+   IF ( return_value == 0 )  THEN
+
+      WRITE(temp_string,*) file_id
+      CALL internal_message( 'debug', routine_name // &
+                                      ': close file (file_id=' // TRIM( temp_string ) // ')' )
+
+      CLOSE( file_id, IOSTAT=return_value )
+      IF ( return_value /= 0 )  THEN
+         WRITE(temp_string,*) file_id
+         CALL internal_message( 'error',        &
+                                routine_name // &
+                                ': cannot close file (file_id=' // TRIM( temp_string ) // ')' )
+      ENDIF
+
+   ENDIF
+
+END SUBROUTINE binary_finalize
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Message routine writing debug information into the debug file
+!> or creating the error message string.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE internal_message( level, string )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  level   !< message importance level
+   CHARACTER(LEN=*), INTENT(IN) ::  string  !< message string
+
+
+   IF ( TRIM( level ) == 'error' )  THEN
+
+      WRITE( internal_error_message, '(A,A)' ) 'DOM ERROR: ', string
+
+   ELSEIF ( TRIM( level ) == 'debug'  .AND.  print_debug_output )  THEN
+
+      WRITE( debug_output_unit, '(A,A)' ) 'DOM DEBUG: ', string
+      FLUSH( debug_output_unit )
+
+   ENDIF
+
+END SUBROUTINE internal_message
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Return the last created error message.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE binary_get_error_message( error_message )
+
+   CHARACTER(LEN=800), INTENT(OUT) ::  error_message  !< return error message to main program
+
+
+   error_message = internal_error_message
+
+END SUBROUTINE binary_get_error_message
+
+
+END MODULE data_output_binary_module
Index: palm/trunk/SOURCE/data_output_module.f90
===================================================================
--- palm/trunk/SOURCE/data_output_module.f90	(revision 4070)
+++ palm/trunk/SOURCE/data_output_module.f90	(revision 4070)
@@ -0,0 +1,3506 @@
+!> @file data_output_module.f90
+!--------------------------------------------------------------------------------------------------!
+! This file is part of the PALM model system.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 2019-2019 Leibniz Universitaet Hannover
+!--------------------------------------------------------------------------------------------------!
+!
+! Current revisions:
+! ------------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
+! Initial revision
+!
+!
+! Authors:
+! --------
+!> @author Tobias Gronemeier
+!> @author Helge Knoop
+!
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Data-output module to handle output of variables into output files.
+!>
+!> The module first creates an interal database containing all meta data of all
+!> output quantities. Output files are then inititialized and prepared for
+!> storing data, which are finally written to file.
+!>
+!> @todo Convert variable if type of given values do not fit specified type.
+!> @todo Remove unused variables
+!> @todo How to deal with definition calls after dom_start_output is called? Should it be allowed
+!>       to define new files after that (which is technically possible)?
+!> @todo Remove iwp from index (and similar) variables.
+!--------------------------------------------------------------------------------------------------!
+MODULE data_output_module
+
+   USE kinds
+
+   USE data_output_netcdf4_serial_module, &
+      ONLY: netcdf4_serial_init_dimension, &
+            netcdf4_serial_get_error_message, &
+            netcdf4_serial_init_end, &
+            netcdf4_serial_init_module, &
+            netcdf4_serial_init_variable, &
+            netcdf4_serial_finalize, &
+            netcdf4_serial_open_file, &
+            netcdf4_serial_write_attribute, &
+            netcdf4_serial_write_variable
+
+   USE data_output_netcdf4_parallel_module, &
+      ONLY: netcdf4_parallel_init_dimension, &
+            netcdf4_parallel_get_error_message, &
+            netcdf4_parallel_init_end, &
+            netcdf4_parallel_init_module, &
+            netcdf4_parallel_init_variable, &
+            netcdf4_parallel_finalize, &
+            netcdf4_parallel_open_file, &
+            netcdf4_parallel_write_attribute, &
+            netcdf4_parallel_write_variable
+
+   USE data_output_binary_module, &
+      ONLY: binary_finalize, &
+            binary_get_error_message, &
+            binary_init_dimension, &
+            binary_init_end, &
+            binary_init_module, &
+            binary_init_variable, &
+            binary_open_file, &
+            binary_write_attribute, &
+            binary_write_variable
+
+   IMPLICIT NONE
+
+   INTEGER(iwp), PARAMETER ::  charlen = 100_iwp  !< maximum length of character variables
+
+   TYPE attribute_type
+      CHARACTER(LEN=charlen) ::  data_type = ''  !< data type
+      CHARACTER(LEN=charlen) ::  name            !< attribute name
+      CHARACTER(LEN=charlen) ::  value_char      !< attribute value if character
+      INTEGER(KIND=1)        ::  value_int8      !< attribute value if 8bit integer
+      INTEGER(KIND=2)        ::  value_int16     !< attribute value if 16bit integer
+      INTEGER(KIND=4)        ::  value_int32     !< attribute value if 32bit integer
+      REAL(KIND=4)           ::  value_real32    !< attribute value if 32bit real
+      REAL(KIND=8)           ::  value_real64    !< attribute value if 64bit real
+   END TYPE attribute_type
+
+   TYPE variable_type
+      CHARACTER(LEN=charlen) ::  data_type = ''       !< data type
+      CHARACTER(LEN=charlen) ::  name                 !< variable name
+      INTEGER(iwp)           ::  id = 0               !< id within file
+      LOGICAL                ::  is_global = .FALSE.  !< true if global variable
+      LOGICAL                ::  is_init = .FALSE.    !< true if initialized
+      CHARACTER(LEN=charlen), DIMENSION(:), ALLOCATABLE ::  dimension_names  !< list of dimension names
+      INTEGER(iwp),           DIMENSION(:), ALLOCATABLE ::  dimension_ids    !< list of dimension ids
+      TYPE(attribute_type),   DIMENSION(:), ALLOCATABLE ::  attributes       !< list of attributes
+   END TYPE variable_type
+
+   TYPE dimension_type
+      CHARACTER(LEN=charlen) ::  data_type = ''       !< data type
+      CHARACTER(LEN=charlen) ::  name                 !< dimension name
+      INTEGER(iwp)           ::  id = 0               !< dimension id within file
+      INTEGER(iwp)           ::  length               !< length of dimension
+      INTEGER(iwp)           ::  length_mask          !< length of masked dimension
+      INTEGER(iwp)           ::  var_id = 0           !< associated variable id within file
+      LOGICAL                ::  is_init = .FALSE.    !< true if initialized
+      LOGICAL                ::  is_masked = .FALSE.  !< true if masked
+      INTEGER(iwp),    DIMENSION(2)              ::  bounds                !< lower and upper bound of dimension
+      INTEGER(iwp),    DIMENSION(:), ALLOCATABLE ::  masked_indices        !< list of masked indices of dimension
+      INTEGER(KIND=1), DIMENSION(:), ALLOCATABLE ::  masked_values_int8    !< masked dimension values if 16bit integer
+      INTEGER(KIND=2), DIMENSION(:), ALLOCATABLE ::  masked_values_int16   !< masked dimension values if 16bit integer
+      INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE ::  masked_values_int32   !< masked dimension values if 32bit integer
+      INTEGER(iwp),    DIMENSION(:), ALLOCATABLE ::  masked_values_intwp   !< masked dimension values if working-precision integer
+      INTEGER(KIND=1), DIMENSION(:), ALLOCATABLE ::  values_int8           !< dimension values if 16bit integer
+      INTEGER(KIND=2), DIMENSION(:), ALLOCATABLE ::  values_int16          !< dimension values if 16bit integer
+      INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE ::  values_int32          !< dimension values if 32bit integer
+      INTEGER(iwp),    DIMENSION(:), ALLOCATABLE ::  values_intwp          !< dimension values if working-precision integer
+      LOGICAL,         DIMENSION(:), ALLOCATABLE ::  mask                  !< mask
+      REAL(KIND=4),    DIMENSION(:), ALLOCATABLE ::  masked_values_real32  !< masked dimension values if 32bit real
+      REAL(KIND=8),    DIMENSION(:), ALLOCATABLE ::  masked_values_real64  !< masked dimension values if 64bit real
+      REAL(wp),        DIMENSION(:), ALLOCATABLE ::  masked_values_realwp  !< masked dimension values if working-precision real
+      REAL(KIND=4),    DIMENSION(:), ALLOCATABLE ::  values_real32         !< dimension values if 32bit real
+      REAL(KIND=8),    DIMENSION(:), ALLOCATABLE ::  values_real64         !< dimension values if 64bit real
+      REAL(wp),        DIMENSION(:), ALLOCATABLE ::  values_realwp         !< dimension values if working-precision real
+      TYPE(attribute_type), DIMENSION(:), ALLOCATABLE ::  attributes       !< list of attributes
+   END TYPE dimension_type
+
+   TYPE file_type
+      CHARACTER(LEN=charlen) ::  format = ''        !< file format
+      CHARACTER(LEN=charlen) ::  name = ''          !< file name
+      INTEGER(iwp)           ::  id = 0             !< id of file
+      LOGICAL                ::  is_init = .FALSE.  !< true if initialized
+      TYPE(attribute_type), DIMENSION(:), ALLOCATABLE ::  attributes  !< list of attributes
+      TYPE(dimension_type), DIMENSION(:), ALLOCATABLE ::  dimensions  !< list of dimensions
+      TYPE(variable_type),  DIMENSION(:), ALLOCATABLE ::  variables   !< list of variables
+   END TYPE file_type
+
+
+   CHARACTER(LEN=charlen) ::  output_file_format = 'binary'  !< file format (namelist parameter)
+
+   CHARACTER(LEN=800) ::  internal_error_message = '' !< string containing the last error message
+   CHARACTER(LEN=800) ::  temp_string                 !< dummy string
+
+   INTEGER(iwp) ::  debug_output_unit  !< Fortran Unit Number of the debug-output file
+   INTEGER(iwp) ::  nf                 !< number of files
+
+   INTEGER(iwp), PARAMETER ::  no_var_id = -1  !< value of var_id if no variable is selected
+
+   LOGICAL ::  print_debug_output = .FALSE.  !< if true, debug output is printed
+
+   TYPE(file_type), DIMENSION(:), ALLOCATABLE ::  files  !< file list
+
+   SAVE
+
+   PRIVATE
+
+   !> Initialize the data-output module
+   INTERFACE dom_init
+      MODULE PROCEDURE dom_init
+   END INTERFACE dom_init
+
+   !> Add files to database
+   INTERFACE dom_def_file
+      MODULE PROCEDURE dom_def_file
+   END INTERFACE dom_def_file
+
+   !> Add dimensions to database
+   INTERFACE dom_def_dim
+      MODULE PROCEDURE dom_def_dim
+   END INTERFACE dom_def_dim
+
+   !> Add variables to database
+   INTERFACE dom_def_var
+      MODULE PROCEDURE dom_def_var
+   END INTERFACE dom_def_var
+
+   !> Add attributes to database
+   INTERFACE dom_def_att
+      MODULE PROCEDURE dom_def_att_char
+      MODULE PROCEDURE dom_def_att_int8
+      MODULE PROCEDURE dom_def_att_int16
+      MODULE PROCEDURE dom_def_att_int32
+      MODULE PROCEDURE dom_def_att_real32
+      MODULE PROCEDURE dom_def_att_real64
+   END INTERFACE dom_def_att
+
+   !> Prepare for output: evaluate database and create files
+   INTERFACE dom_start_output
+      MODULE PROCEDURE dom_start_output
+   END INTERFACE dom_start_output
+
+   !> Write variables to file
+   INTERFACE dom_write_var
+      MODULE PROCEDURE dom_write_var
+   END INTERFACE dom_write_var
+
+   !> Last actions required for output befor termination
+   INTERFACE dom_finalize_output
+      MODULE PROCEDURE dom_finalize_output
+   END INTERFACE dom_finalize_output
+
+   !> Return error message
+   INTERFACE dom_get_error_message
+      MODULE PROCEDURE dom_get_error_message
+   END INTERFACE dom_get_error_message
+
+   PUBLIC &
+      dom_database_debug_output, &
+      dom_def_att, &
+      dom_def_dim, &
+      dom_def_file, &
+      dom_def_var, &
+      dom_finalize_output, &
+      dom_get_error_message, &
+      dom_init, &
+      dom_start_output, &
+      dom_write_var
+
+CONTAINS
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize data-output module
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE dom_init( program_debug_output_unit, debug_output )
+
+   INTEGER(iwp), INTENT(IN) ::  program_debug_output_unit  !< file unit number for debug output
+
+   LOGICAL, INTENT(IN) ::  debug_output  !< if true, debug output is printed
+
+
+   debug_output_unit = program_debug_output_unit
+
+   print_debug_output = debug_output
+
+   CALL binary_init_module( debug_output_unit, debug_output, no_var_id )
+
+   CALL netcdf4_serial_init_module( debug_output_unit, debug_output, no_var_id )
+
+   CALL netcdf4_parallel_init_module( debug_output_unit, debug_output, no_var_id )
+
+END SUBROUTINE dom_init
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Debugging output. Print contents of output database to terminal.
+!>
+!> @bug Routine currently not working properly due to changes in dimension/attribute value types.
+!> @todo Call routine via debug flag.
+!> @todo Use dim and var within loops to enhance readability of the code.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE dom_database_debug_output
+
+   ! CHARACTER(LEN=100) ::  dim_string  !< list of dimension names as single string
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_database_debug_output'  !< name of this routine
+
+   ! INTEGER(iwp) ::  d     !< loop index
+   ! INTEGER(iwp) ::  f     !< loop index
+   ! INTEGER(iwp) ::  i     !< loop index
+   ! INTEGER(iwp) ::  ndim  !< number of dimensions
+   !
+   ! TYPE(dimension_type) ::  dim  !< dimension in file
+   !
+   ! TYPE(variable_type) ::  var   !< variable in file
+
+
+   CALL internal_message( 'error', routine_name // ': routine currently not available' )
+
+   ! nf = SIZE(files)
+   !
+   ! WRITE(*, '(A,I3)') '    number of files: ', nf
+   ! WRITE(*, '(A)')    '    files:'
+   ! DO  f = 1, nf
+   !    WRITE(*,'(//6X,A)') 'name: '//TRIM(files(f)%name)
+   !    WRITE(*,'(6X,A)') 'attributes:'
+   !    IF ( ALLOCATED(files(f)%attributes) )  THEN
+   !       CALL print_attributes( files(f)%attributes )
+   !    ELSE
+   !       WRITE(*,'(8X,A)') '--none--'
+   !    ENDIF
+   !
+   !    WRITE(*,'(/6X,A)') 'dimensions:'
+   !    IF ( ALLOCATED(files(f)%dimensions )  THEN
+   !       ndim = SIZE(files(f)%dimensions)
+   !       DO  d = 1, ndim
+   !
+   !          WRITE(*,'(/8X,A,I4/8X,A)') &
+   !             'name: '//TRIM(files(f)%dimensions(d)%name)//',  '// &
+   !             'length: ', files(f)%dimensions(d)%length, &
+   !             'values:'
+   !          IF ( ALLOCATED(files(f)%dimensions(d)%values_int32) )  THEN
+   !             WRITE(*,'(10X,5(I5,1X),A)') files(f)%dimensions(d)%values_int32(0:MIN(4,files(f)%dimensions(d)%length)), '...'
+   !          ELSE
+   !             WRITE(*,'(10X,5(F8.2,1X),A)') files(f)%dimensions(d)%values_real64(0:MIN(4,files(f)%dimensions(d)%length)), '...'
+   !          ENDIF
+   !          IF ( ALLOCATED(files(f)%dimensions(d)%mask) )  THEN
+   !             WRITE(*,'(8X,A)') 'mask:'
+   !             WRITE(*,'(10X,5(L,1X),A)') files(f)%dimensions(d)%mask(0:MIN(4,files(f)%dimensions(d)%length)), '...'
+   !          ENDIF
+   !          WRITE(*,'(8X,A)') 'attributes:'
+   !          IF ( ALLOCATED( files(f)%dimensions(d)%attributes ) )  THEN
+   !             CALL print_attributes( files(f)%dimensions(d)%attributes )
+   !          ELSE
+   !             WRITE(*,'(10X,A)') '--none--'
+   !          ENDIF
+   !
+   !       ENDDO
+   !    ELSE
+   !       WRITE(*,'(8X,A)') '--none--'
+   !    ENDIF
+   !
+   !    WRITE(*,'(/6X,A)') 'variables:'
+   !    IF ( ALLOCATED(files(f)%variables) )  THEN
+   !       ndim = SIZE(files(f)%variables)
+   !       DO  d = 1, ndim
+   !          dim_string = '('
+   !          DO  i = 1, SIZE(files(f)%variables(d)%dimension_names)
+   !             dim_string = TRIM(dim_string)//' '//TRIM(files(f)%variables(d)%dimension_names(i))
+   !          ENDDO
+   !          dim_string = TRIM(dim_string)//')'
+   !          WRITE(*,'(/8X,A)') &
+   !             'name: '//TRIM(files(f)%variables(d)%name)//TRIM(dim_string)
+   !          WRITE(*,'(/8X,A)') &
+   !             'type: '//TRIM(files(f)%variables(d)%data_type)
+   !          WRITE(*,'(8X,A)') &
+   !             'ID: '
+   !          WRITE(*,'(8X,A)') 'attributes:'
+   !          IF ( ALLOCATED( files(f)%variables(d)%attributes ) )  THEN
+   !             CALL print_attributes( files(f)%variables(d)%attributes )
+   !          ELSE
+   !             WRITE(*,'(10X,A)') '--none--'
+   !          ENDIF
+   !       ENDDO
+   !    ELSE
+   !       WRITE(*,'(8X,A)') '--none--'
+   !    ENDIF
+   !
+   !    WRITE(*,'(4X,40("#"))')
+   !
+   ! ENDDO
+   !
+   ! CONTAINS
+   !
+   !    !------------------------------------------------------------------------!
+   !    ! Description:
+   !    ! ------------
+   !    !> Print given list of attributes.
+   !    !------------------------------------------------------------------------!
+   !    SUBROUTINE print_attributes( attributes )
+   !
+   !       INTEGER(iwp) ::  a     !< loop index
+   !       INTEGER(iwp) ::  natt  !< number of attributes
+   !
+   !       TYPE(attribute_type), DIMENSION(:), INTENT(IN) ::  attributes  !< list of attributes
+   !
+   !       natt = SIZE(attributes)
+   !       DO  a = 1, natt
+   !          IF ( TRIM(attributes(a)%data_type) == 'int32' )  THEN
+   !             WRITE(*,'(10X,A,1X,I10)') &
+   !                '"'//TRIM(attributes(a)%name)//'" = ', &
+   !                attributes(a)%value_int32
+   !          ELSEIF (  TRIM(attributes(a)%data_type) == 'real32' )  THEN
+   !             WRITE(*,'(10X,A,1X,F10.3)') &
+   !                '"'//TRIM(attributes(a)%name)//'" = ', &
+   !                attributes(a)%value_real64
+   !          ELSEIF (  TRIM(attributes(a)%data_type) == 'char' )  THEN
+   !             WRITE(*,'(10X,A,1X,A)') &
+   !                '"'//TRIM(attributes(a)%name)//'" = ', &
+   !                '"'//TRIM(attributes(a)%value_char)//'"'
+   !          ENDIF
+   !       ENDDO
+   !
+   !    END SUBROUTINE print_attributes
+
+END SUBROUTINE dom_database_debug_output
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Define output file.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION dom_def_file( filename, format ) RESULT( return_value )
+
+   CHARACTER(LEN=*), INTENT(IN)           ::  filename  !< name of file to be created
+   CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  format    !< format of file to be created
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_def_file'  !< name of this routine
+
+   INTEGER(iwp) ::  f             !< loop index
+   INTEGER(iwp) ::  return_value  !< return value
+
+   TYPE(file_type), DIMENSION(:), ALLOCATABLE ::  files_tmp  !< temporary file list
+
+
+   return_value = 0
+
+   !-- Allocate file list or extend it by 1
+   IF ( .NOT. ALLOCATED( files ) ) THEN
+
+      nf = 1
+      ALLOCATE( files(nf) )
+
+   ELSE
+
+      nf = SIZE( files )
+      !-- Check if file already exists
+      DO  f = 1, nf
+         IF ( files(f)%name == TRIM( filename ) )  THEN
+            return_value = 1
+            CALL internal_message( 'error', routine_name // ': file "' // TRIM( filename ) // &
+                                            '" already exists' )
+            EXIT
+         ENDIF
+      ENDDO
+
+      !-- Extend file list
+      IF ( return_value == 0 )  THEN
+         ALLOCATE( files_tmp(nf) )
+         files_tmp = files
+         DEALLOCATE( files )
+         nf = nf + 1
+         ALLOCATE( files(nf) )
+         files(:nf-1) = files_tmp
+         DEALLOCATE( files_tmp )
+      ENDIF
+
+   ENDIF
+
+   !-- Add new file to database
+   IF ( return_value == 0 )  THEN
+      files(nf)%name = TRIM( filename )
+      IF ( PRESENT( format ) )  THEN
+         files(nf)%format = TRIM( format )
+      ELSE
+         files(nf)%format = TRIM( output_file_format )
+      ENDIF
+   ENDIF
+
+END FUNCTION dom_def_file
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Define dimension of type integer.
+!>
+!> @todo Convert given values into selected output_type.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION dom_def_dim( filename, name, output_type, bounds,                   &
+                      values_int8, values_int16, values_int32, values_intwp, &
+                      values_real32, values_real64, values_realwp,           &
+                      mask ) RESULT( return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  filename     !< name of file
+   CHARACTER(LEN=*), INTENT(IN) ::  name         !< name of dimension
+   CHARACTER(LEN=*), INTENT(IN) ::  output_type  !< data type of dimension variable in output file
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_def_dim'  !< name of this routine
+
+   INTEGER(iwp) ::  d             !< loop index
+   INTEGER(iwp) ::  f             !< loop index
+   INTEGER(iwp) ::  i             !< loop index
+   INTEGER(iwp) ::  j             !< loop index
+   INTEGER(iwp) ::  ndim          !< number of dimensions in file
+   INTEGER(iwp) ::  return_value  !< return value
+
+   INTEGER(iwp), DIMENSION(:), INTENT(IN) ::  bounds  !< lower and upper bound of dimension variable
+
+   INTEGER(KIND=1), DIMENSION(:), INTENT(IN), OPTIONAL ::  values_int8   !< values of dimension
+   INTEGER(KIND=2), DIMENSION(:), INTENT(IN), OPTIONAL ::  values_int16  !< values of dimension
+   INTEGER(KIND=4), DIMENSION(:), INTENT(IN), OPTIONAL ::  values_int32  !< values of dimension
+   INTEGER(iwp),    DIMENSION(:), INTENT(IN), OPTIONAL ::  values_intwp  !< values of dimension
+
+   LOGICAL, DIMENSION(:), INTENT(IN), OPTIONAL ::  mask  !< mask of dimesion
+
+   REAL(KIND=4), DIMENSION(:), INTENT(IN), OPTIONAL ::  values_real32  !< values of dimension
+   REAL(KIND=8), DIMENSION(:), INTENT(IN), OPTIONAL ::  values_real64  !< values of dimension
+   REAL(wp),     DIMENSION(:), INTENT(IN), OPTIONAL ::  values_realwp  !< values of dimension
+
+   TYPE(dimension_type) ::  dimension  !< new dimension
+
+   TYPE(dimension_type), DIMENSION(:), ALLOCATABLE ::  dims_tmp  !< temporary dimension list
+
+
+   return_value = 0
+
+   dimension%name      = TRIM( name )
+   dimension%data_type = TRIM( output_type )
+
+   !-- Check dimension bounds and allocate dimension according to bounds
+   IF ( SIZE( bounds ) == 1 )  THEN
+
+      !-- Dimension has only lower bound, which means it changes its size
+      !-- during simulation.
+      !-- Set length to -1 as indicator.
+      dimension%bounds(:) = bounds(1)
+      dimension%length    = -1_iwp
+
+      IF ( PRESENT( mask ) )  THEN
+         return_value = 1
+         CALL internal_message( 'error', routine_name //                            &
+                                         ': unlimited dimension' // TRIM( name ) // &
+                                         ' in file ' // TRIM( filename ) // ' cannot be masked' )
+      ENDIF
+
+   ELSEIF ( SIZE( bounds ) == 2 )  THEN
+
+      dimension%bounds = bounds
+      dimension%length = bounds(2) - bounds(1) + 1
+
+      !-- Save dimension values
+      IF ( PRESENT( values_int8 ) )  THEN
+         ALLOCATE( dimension%values_int8(dimension%bounds(1):dimension%bounds(2)) )
+         IF ( SIZE( values_int8 ) == dimension%length )  THEN
+            dimension%values_int8 = values_int8
+         ELSEIF ( SIZE( values_int8 ) == 1 )  THEN
+            dimension%values_int8(:) = values_int8
+         ELSE
+            return_value = 2
+         ENDIF
+      ELSEIF( PRESENT( values_int16 ) )  THEN
+         ALLOCATE( dimension%values_int16(dimension%bounds(1):dimension%bounds(2)) )
+         IF ( SIZE( values_int16 ) == dimension%length )  THEN
+            dimension%values_int16 = values_int16
+         ELSEIF ( SIZE( values_int16 ) == 1 )  THEN
+            dimension%values_int16(:) = values_int16
+         ELSE
+            return_value = 2
+         ENDIF
+      ELSEIF( PRESENT( values_int32 ) )  THEN
+         ALLOCATE( dimension%values_int32(dimension%bounds(1):dimension%bounds(2)) )
+         IF ( SIZE( values_int32 ) == dimension%length )  THEN
+            dimension%values_int32 = values_int32
+         ELSEIF ( SIZE( values_int32 ) == 1 )  THEN
+            dimension%values_int32(:) = values_int32
+         ELSE
+            return_value = 2
+         ENDIF
+      ELSEIF( PRESENT( values_intwp ) )  THEN
+         ALLOCATE( dimension%values_intwp(dimension%bounds(1):dimension%bounds(2)) )
+         IF ( SIZE( values_intwp ) == dimension%length )  THEN
+            dimension%values_intwp = values_intwp
+         ELSEIF ( SIZE( values_intwp ) == 1 )  THEN
+            dimension%values_intwp(:) = values_intwp
+         ELSE
+            return_value = 2
+         ENDIF
+      ELSEIF( PRESENT( values_real32 ) )  THEN
+         ALLOCATE( dimension%values_real32(dimension%bounds(1):dimension%bounds(2)) )
+         IF ( SIZE( values_real32 ) == dimension%length )  THEN
+            dimension%values_real32 = values_real32
+         ELSEIF ( SIZE( values_real32 ) == 1 )  THEN
+            dimension%values_real32(:) = values_real32
+         ELSE
+            return_value = 2
+         ENDIF
+      ELSEIF( PRESENT( values_real64 ) )  THEN
+         ALLOCATE( dimension%values_real64(dimension%bounds(1):dimension%bounds(2)) )
+         IF ( SIZE( values_real64 ) == dimension%length )  THEN
+            dimension%values_real64 = values_real64
+         ELSEIF ( SIZE( values_real64 ) == 1 )  THEN
+            dimension%values_real64(:) = values_real64
+         ELSE
+            return_value = 2
+         ENDIF
+      ELSEIF( PRESENT( values_realwp ) )  THEN
+         ALLOCATE( dimension%values_realwp(dimension%bounds(1):dimension%bounds(2)) )
+         IF ( SIZE( values_realwp ) == dimension%length )  THEN
+            dimension%values_realwp = values_realwp
+         ELSEIF ( SIZE( values_realwp ) == 1 )  THEN
+            dimension%values_realwp(:) = values_realwp
+         ELSE
+            return_value = 2
+         ENDIF
+      ELSE
+         return_value = 1
+         CALL internal_message( 'error', routine_name // ': ' // &
+                                         TRIM( name ) // ': no values given' )
+      ENDIF
+
+      IF ( return_value == 2 )  THEN
+         return_value = 1
+         CALL internal_message( 'error', routine_name //                   &
+                                         ': dimension ' // TRIM( name ) // &
+                                         ': number of values and given bounds do not match' )
+      ENDIF
+
+      !-- Initialize mask
+      IF ( PRESENT( mask )  .AND.  return_value == 0 )  THEN
+
+         dimension%is_masked = .TRUE.
+
+         IF ( dimension%length == SIZE( mask ) )  THEN
+
+            dimension%length_mask = COUNT( mask )
+
+            ALLOCATE( dimension%mask(dimension%bounds(1):dimension%bounds(2)) )
+            ALLOCATE( dimension%masked_indices(0:dimension%length_mask-1) )
+
+            dimension%mask = mask
+
+            !-- Save masked positions and masked values
+            IF ( ALLOCATED( dimension%values_int8 ) )  THEN
+
+               ALLOCATE( dimension%masked_values_int8(0:dimension%length_mask-1) )
+               j = 0
+               DO  i = 0, dimension%length-1
+                  IF ( dimension%mask(i) )  THEN
+                     dimension%masked_values_int8(j) = dimension%values_int8(i)
+                     dimension%masked_indices(j) = i
+                     j = j + 1
+                  ENDIF
+               ENDDO
+
+            ELSEIF ( ALLOCATED( dimension%values_int16 ) )  THEN
+
+               ALLOCATE( dimension%masked_values_int16(0:dimension%length_mask-1) )
+               j = 0
+               DO  i = 0, dimension%length-1
+                  IF ( dimension%mask(i) )  THEN
+                     dimension%masked_values_int16(j) = dimension%values_int16(i)
+                     dimension%masked_indices(j) = i
+                     j = j + 1
+                  ENDIF
+               ENDDO
+
+            ELSEIF ( ALLOCATED( dimension%values_int32 ) )  THEN
+
+               ALLOCATE( dimension%masked_values_int32(0:dimension%length_mask-1) )
+               j = 0
+               DO  i = 0, dimension%length-1
+                  IF ( dimension%mask(i) )  THEN
+                     dimension%masked_values_int32(j) = dimension%values_int32(i)
+                     dimension%masked_indices(j) = i
+                     j = j + 1
+                  ENDIF
+               ENDDO
+
+            ELSEIF ( ALLOCATED( dimension%values_intwp ) )  THEN
+
+               ALLOCATE( dimension%masked_values_intwp(0:dimension%length_mask-1) )
+               j = 0
+               DO  i = 0, dimension%length-1
+                  IF ( dimension%mask(i) )  THEN
+                     dimension%masked_values_intwp(j) = dimension%values_intwp(i)
+                     dimension%masked_indices(j) = i
+                     j = j + 1
+                  ENDIF
+               ENDDO
+
+            ELSEIF ( ALLOCATED( dimension%values_real32 ) )  THEN
+
+               ALLOCATE( dimension%masked_values_real32(0:dimension%length_mask-1) )
+               j = 0
+               DO  i = 0, dimension%length-1
+                  IF ( dimension%mask(i) )  THEN
+                     dimension%masked_values_real32(j) = dimension%values_real32(i)
+                     dimension%masked_indices(j) = i
+                     j = j + 1
+                  ENDIF
+               ENDDO
+
+            ELSEIF ( ALLOCATED(dimension%values_real64) )  THEN
+
+               ALLOCATE( dimension%masked_values_real64(0:dimension%length_mask-1) )
+               j = 0
+               DO  i = 0, dimension%length-1
+                  IF ( dimension%mask(i) )  THEN
+                     dimension%masked_values_real64(j) = dimension%values_real64(i)
+                     dimension%masked_indices(j) = i
+                     j = j + 1
+                  ENDIF
+               ENDDO
+
+            ELSEIF ( ALLOCATED(dimension%values_realwp) )  THEN
+
+               ALLOCATE( dimension%masked_values_realwp(0:dimension%length_mask-1) )
+               j = 0
+               DO  i = 0, dimension%length-1
+                  IF ( dimension%mask(i) )  THEN
+                     dimension%masked_values_realwp(j) = dimension%values_realwp(i)
+                     dimension%masked_indices(j) = i
+                     j = j + 1
+                  ENDIF
+               ENDDO
+
+            ENDIF
+
+         ELSE
+            return_value = 1
+            CALL internal_message( 'error', routine_name //                   &
+                                            ': dimension ' // TRIM( name ) // &
+                                            ': size of mask and given bounds do not match' )
+         ENDIF
+
+      ENDIF
+
+   ELSE
+
+      return_value = 1
+      CALL internal_message( 'error', routine_name //                   &
+                                      ': dimension ' // TRIM( name ) // &
+                                      ': At least one but no more than two bounds must be given' )
+
+   ENDIF
+
+   !-- Add dimension to database
+   IF ( return_value == 0 )  THEN
+
+      DO  f = 1, nf
+
+         IF ( TRIM( filename ) == files(f)%name )  THEN
+
+            IF ( .NOT. ALLOCATED( files(f)%dimensions ) )  THEN
+
+               ndim = 1
+               ALLOCATE( files(f)%dimensions(ndim) )
+
+            ELSE
+
+               ndim = SIZE( files(f)%dimensions )
+
+               !-- Check if dimension already exists in file
+               DO  d = 1, ndim
+                  IF ( files(f)%dimensions(d)%name == dimension%name )  THEN
+                     return_value = 1
+                     CALL internal_message( 'error',            &
+                             routine_name //                    &
+                             ': dimension "' // TRIM( name ) // &
+                             '" already exists in file "' // TRIM( filename ) // '"' )
+                     EXIT
+                  ENDIF
+               ENDDO
+
+               !-- Extend dimension list
+               IF ( return_value == 0 )  THEN
+                  ALLOCATE( dims_tmp(ndim) )
+                  dims_tmp = files(f)%dimensions
+                  DEALLOCATE( files(f)%dimensions )
+                  ndim = ndim + 1
+                  ALLOCATE( files(f)%dimensions(ndim) )
+                  files(f)%dimensions(:ndim-1) = dims_tmp
+                  DEALLOCATE( dims_tmp )
+               ENDIF
+
+            ENDIF
+
+            !-- Add new dimension to database
+            IF ( return_value == 0 )  files(f)%dimensions(ndim) = dimension
+
+            EXIT
+
+         ENDIF
+      ENDDO
+
+      IF ( f > nf )  THEN
+         return_value = 1
+         CALL internal_message( 'error', routine_name // &
+                                         ': file "' // TRIM( filename ) // '" not found' )
+      ENDIF
+
+   ENDIF
+
+END FUNCTION dom_def_dim
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Add variable to database.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION dom_def_var( filename, name, dimension_names, output_type, is_global ) &
+            RESULT( return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  filename     !< name of file
+   CHARACTER(LEN=*), INTENT(IN) ::  name         !< name of variable
+   CHARACTER(LEN=*), INTENT(IN) ::  output_type  !< data type of variable
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_def_var'  !< name of this routine
+
+   CHARACTER(LEN=*), DIMENSION(:), INTENT(IN) ::  dimension_names  !< list of dimension names
+
+   INTEGER(iwp) ::  d             !< loop index
+   INTEGER(iwp) ::  f             !< loop index
+   INTEGER(iwp) ::  i             !< loop index
+   INTEGER(iwp) ::  nvar          !< number of variables in file
+   INTEGER(iwp) ::  return_value  !< return value
+
+   LOGICAL                       ::  found = .FALSE.  !< true if requested dimension is defined in file
+   LOGICAL, INTENT(IN), OPTIONAL ::  is_global        !< true if variable is global (same on all PE)
+
+   TYPE(variable_type) ::  variable  !< new variable
+
+   TYPE(variable_type), DIMENSION(:), ALLOCATABLE ::  vars_tmp  !< temporary variable list
+
+
+   return_value = 0
+
+   variable%name = TRIM( name )
+
+   ALLOCATE( variable%dimension_names(SIZE( dimension_names )) )
+   ALLOCATE( variable%dimension_ids(SIZE( dimension_names )) )
+
+   variable%dimension_names = dimension_names
+   variable%data_type = TRIM( output_type )
+
+   IF ( PRESENT( is_global ) )  THEN
+      variable%is_global = is_global
+   ELSE
+      variable%is_global = .FALSE.
+   ENDIF
+
+   !-- Add variable to database
+   DO  f = 1, nf
+
+      IF ( TRIM( filename ) == files(f)%name )  THEN
+
+         !-- Check if dimensions assigned to variable are defined within file
+         IF ( ALLOCATED( files(f)%dimensions ) )  THEN
+
+            DO  i = 1, SIZE( variable%dimension_names )
+               found = .FALSE.
+               DO  d = 1, SIZE( files(f)%dimensions )
+                  IF ( files(f)%dimensions(d)%name == variable%dimension_names(i) )  THEN
+                     found = .TRUE.
+                     EXIT
+                  ENDIF
+               ENDDO
+               IF ( .NOT. found )  THEN
+                  return_value = 1
+                  CALL internal_message( 'error', &
+                                         routine_name //                      &
+                                         ': variable "' // TRIM( name ) //    &
+                                         '" in file "' // TRIM( filename ) // &
+                                         '": required dimension "' //         &
+                                         TRIM( variable%dimension_names(i) ) // '" not defined' )
+                  EXIT
+               ENDIF
+            ENDDO
+
+         ELSE
+
+            return_value = 1
+            CALL internal_message( 'error', routine_name //                   &
+                                            ': file "' // TRIM( filename ) // &
+                                            '" has no dimensions defined' )
+
+         ENDIF
+
+         IF ( return_value == 0 )  THEN
+
+            !-- Check if variable already exists
+            IF ( .NOT. ALLOCATED( files(f)%variables ) )  THEN
+
+               nvar = 1
+               ALLOCATE( files(f)%variables(nvar) )
+
+            ELSE
+
+               nvar = SIZE( files(f)%variables )
+               DO  i = 1, nvar
+                  IF ( files(f)%variables(i)%name == variable%name )  THEN
+                     return_value = 1
+                     CALL internal_message( 'error', routine_name //                      &
+                                                     ': variable "' // TRIM( name ) //    &
+                                                     '" in file "' // TRIM( filename ) // &
+                                                     '": variable already exists' )
+                     EXIT
+                  ENDIF
+               ENDDO
+
+               IF ( return_value == 0 )  THEN
+                  !-- Extend variable list
+                  ALLOCATE( vars_tmp(nvar) )
+                  vars_tmp = files(f)%variables
+                  DEALLOCATE( files(f)%variables )
+                  nvar = nvar + 1
+                  ALLOCATE( files(f)%variables(nvar) )
+                  files(f)%variables(:nvar-1) = vars_tmp
+                  DEALLOCATE( vars_tmp )
+               ENDIF
+
+            ENDIF
+
+            !-- Add new variable to database
+            IF ( return_value == 0 )  files(f)%variables(nvar) = variable
+
+         ENDIF
+
+         EXIT
+
+      ENDIF
+
+   ENDDO
+
+   IF ( f > nf )  THEN
+      return_value = 1
+      CALL internal_message( 'error', routine_name //                    &
+                                      ': variable "' // TRIM( name ) //  &
+                                      '": file "' // TRIM( filename ) // &
+                                      '" not found' )
+   ENDIF
+
+END FUNCTION dom_def_var
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Create attribute with value of type character.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION dom_def_att_char( filename, variable, name, value, append ) RESULT( return_value )
+
+   CHARACTER(LEN=*), INTENT(IN)           ::  filename  !< name of file
+   CHARACTER(LEN=*), INTENT(IN)           ::  name      !< name of attribute
+   CHARACTER(LEN=*), INTENT(IN)           ::  value     !< attribute value
+   CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  variable  !< name of variable
+
+   ! CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_def_att_char'  !< name of routine
+
+   INTEGER(iwp) ::  return_value  !< return value
+
+   LOGICAL                       ::  append_internal  !< same as 'append'
+   LOGICAL, INTENT(IN), OPTIONAL ::  append           !< if true, append value to existing value
+
+   TYPE(attribute_type) ::  attribute  !< new attribute
+
+
+   return_value = 0
+
+   IF ( PRESENT( append ) )  THEN
+      append_internal = append
+   ELSE
+      append_internal = .FALSE.
+   ENDIF
+
+   attribute%name       = TRIM( name )
+   attribute%data_type  = 'char'
+   attribute%value_char = TRIM( value )
+
+   IF ( PRESENT( variable ) )  THEN
+      return_value = dom_def_att_save( TRIM( filename ), TRIM( variable ), &
+                                           attribute=attribute, append=append_internal )
+   ELSE
+      return_value = dom_def_att_save( TRIM( filename ), &
+                                           attribute=attribute, append=append_internal )
+   ENDIF
+
+END FUNCTION dom_def_att_char
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Create attribute with value of type int8.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION dom_def_att_int8( filename, variable, name, value, append ) RESULT( return_value )
+
+   CHARACTER(LEN=*), INTENT(IN)           ::  filename  !< name of file
+   CHARACTER(LEN=*), INTENT(IN)           ::  name      !< name of attribute
+   CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  variable  !< name of variable
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_def_att_int8'  !< name of routine
+
+   INTEGER(KIND=1), INTENT(IN) ::  value  !< attribute value
+
+   INTEGER(iwp) ::  return_value  !< return value
+
+   LOGICAL                       ::  append_internal  !< same as 'append'
+   LOGICAL, INTENT(IN), OPTIONAL ::  append           !< if true, append value to existing value
+
+   TYPE(attribute_type) ::  attribute  !< new attribute
+
+
+   return_value = 0
+
+   IF ( PRESENT( append ) )  THEN
+      IF ( append )  THEN
+         return_value = 1
+         CALL internal_message( 'error',                           &
+                                routine_name //                    &
+                                ': attribute "' // TRIM( name ) // &
+                                '": append of numeric attribute not possible.' )
+      ENDIF
+   ENDIF
+
+   IF ( return_value == 0 )  THEN
+      append_internal = .FALSE.
+
+      attribute%name       = TRIM( name )
+      attribute%data_type  = 'int8'
+      attribute%value_int8 = value
+
+      IF ( PRESENT( variable ) )  THEN
+         return_value = dom_def_att_save( TRIM( filename ), TRIM( variable ), &
+                                              attribute=attribute, append=append_internal )
+      ELSE
+         return_value = dom_def_att_save( TRIM( filename ), &
+                                              attribute=attribute, append=append_internal )
+      ENDIF
+   ENDIF
+
+END FUNCTION dom_def_att_int8
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Create attribute with value of type int16.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION dom_def_att_int16( filename, variable, name, value, append ) RESULT( return_value )
+
+   CHARACTER(LEN=*), INTENT(IN)           ::  filename  !< name of file
+   CHARACTER(LEN=*), INTENT(IN)           ::  name      !< name of attribute
+   CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  variable  !< name of variable
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_def_att_int16'  !< name of routine
+
+   INTEGER(KIND=2), INTENT(IN) ::  value  !< attribute value
+
+   INTEGER(iwp) ::  return_value  !< return value
+
+   LOGICAL                       ::  append_internal  !< same as 'append'
+   LOGICAL, INTENT(IN), OPTIONAL ::  append           !< if true, append value to existing value
+
+   TYPE(attribute_type) ::  attribute  !< new attribute
+
+
+   return_value = 0
+
+   IF ( PRESENT( append ) )  THEN
+      IF ( append )  THEN
+         return_value = 1
+         CALL internal_message( 'error',                           &
+                                routine_name //                    &
+                                ': attribute "' // TRIM( name ) // &
+                                '": append of numeric attribute not possible.' )
+      ENDIF
+   ENDIF
+
+   IF ( return_value == 0 )  THEN
+      append_internal = .FALSE.
+
+      attribute%name        = TRIM( name )
+      attribute%data_type   = 'int16'
+      attribute%value_int16 = value
+
+      IF ( PRESENT( variable ) )  THEN
+         return_value = dom_def_att_save( TRIM( filename ), TRIM( variable ), &
+                                               attribute=attribute, append=append_internal )
+      ELSE
+         return_value = dom_def_att_save( TRIM( filename ), &
+                                               attribute=attribute, append=append_internal )
+      ENDIF
+   ENDIF
+
+END FUNCTION dom_def_att_int16
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Create attribute with value of type int32.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION dom_def_att_int32( filename, variable, name, value, append ) RESULT( return_value )
+
+   CHARACTER(LEN=*), INTENT(IN)           ::  filename  !< name of file
+   CHARACTER(LEN=*), INTENT(IN)           ::  name      !< name of attribute
+   CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  variable  !< name of variable
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_def_att_int32'  !< name of routine
+
+   INTEGER(KIND=4), INTENT(IN) ::  value  !< attribute value
+
+   INTEGER(iwp) ::  return_value  !< return value
+
+   LOGICAL                       ::  append_internal  !< same as 'append'
+   LOGICAL, INTENT(IN), OPTIONAL ::  append           !< if true, append value to existing value
+
+   TYPE(attribute_type) ::  attribute  !< new attribute
+
+
+   return_value = 0
+
+   IF ( PRESENT( append ) )  THEN
+      IF ( append )  THEN
+         return_value = 1
+         CALL internal_message( 'error',                           &
+                                routine_name //                    &
+                                ': attribute "' // TRIM( name ) // &
+                                '": append of numeric attribute not possible.' )
+      ENDIF
+   ENDIF
+
+   IF ( return_value == 0 )  THEN
+      append_internal = .FALSE.
+
+      attribute%name        = TRIM( name )
+      attribute%data_type   = 'int32'
+      attribute%value_int32 = value
+
+      IF ( PRESENT( variable ) )  THEN
+         return_value = dom_def_att_save( TRIM( filename ), TRIM( variable ), &
+                                               attribute=attribute, append=append_internal )
+      ELSE
+         return_value = dom_def_att_save( TRIM( filename ), &
+                                               attribute=attribute, append=append_internal )
+      ENDIF
+   ENDIF
+
+END FUNCTION dom_def_att_int32
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Create attribute with value of type real32.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION dom_def_att_real32( filename, variable, name, value, append ) RESULT( return_value )
+
+   CHARACTER(LEN=*), INTENT(IN)           ::  filename  !< name of file
+   CHARACTER(LEN=*), INTENT(IN)           ::  name      !< name of attribute
+   CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  variable  !< name of variable
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_def_att_real32'  !< name of routine
+
+   INTEGER(iwp) ::  return_value  !< return value
+
+   LOGICAL                       ::  append_internal  !< same as 'append'
+   LOGICAL, INTENT(IN), OPTIONAL ::  append           !< if true, append value to existing value
+
+   REAL(KIND=4), INTENT(IN) ::  value  !< attribute value
+
+   TYPE(attribute_type) ::  attribute  !< new attribute
+
+
+   return_value = 0
+
+   IF ( PRESENT( append ) )  THEN
+      IF ( append )  THEN
+         return_value = 1
+         CALL internal_message( 'error',                           &
+                                routine_name //                    &
+                                ': attribute "' // TRIM( name ) // &
+                                '": append of numeric attribute not possible.' )
+      ENDIF
+   ENDIF
+
+   IF ( return_value == 0 )  THEN
+      append_internal = .FALSE.
+
+      attribute%name         = TRIM( name )
+      attribute%data_type    = 'real32'
+      attribute%value_real32 = value
+
+      IF ( PRESENT( variable ) )  THEN
+         return_value = dom_def_att_save( TRIM( filename ), TRIM( variable ), &
+                                                attribute=attribute, append=append_internal )
+      ELSE
+         return_value = dom_def_att_save( TRIM( filename ), &
+                                                attribute=attribute, append=append_internal )
+      ENDIF
+   ENDIF
+
+END FUNCTION dom_def_att_real32
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Create attribute with value of type real64.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION dom_def_att_real64( filename, variable, name, value, append ) RESULT( return_value )
+
+   CHARACTER(LEN=*), INTENT(IN)           ::  filename  !< name of file
+   CHARACTER(LEN=*), INTENT(IN)           ::  name      !< name of attribute
+   CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  variable  !< name of variable
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_def_att_real64'  !< name of routine
+
+   INTEGER(iwp) ::  return_value  !< return value
+
+   LOGICAL                       ::  append_internal  !< same as 'append'
+   LOGICAL, INTENT(IN), OPTIONAL ::  append           !< if true, append value to existing value
+
+   REAL(KIND=8), INTENT(IN) ::  value  !< attribute value
+
+   TYPE(attribute_type) ::  attribute  !< new attribute
+
+
+   return_value = 0
+
+   IF ( PRESENT( append ) )  THEN
+      IF ( append )  THEN
+         return_value = 1
+         CALL internal_message( 'error',                           &
+                                routine_name //                    &
+                                ': attribute "' // TRIM( name ) // &
+                                '": append of numeric attribute not possible.' )
+      ENDIF
+   ENDIF
+
+   IF ( return_value == 0 )  THEN
+      append_internal = .FALSE.
+
+      attribute%name         = TRIM( name )
+      attribute%data_type    = 'real64'
+      attribute%value_real64 = value
+
+      IF ( PRESENT( variable ) )  THEN
+         return_value = dom_def_att_save( TRIM( filename ), TRIM( variable ), &
+                                                attribute=attribute, append=append_internal )
+      ELSE
+         return_value = dom_def_att_save( TRIM( filename ), &
+                                                attribute=attribute, append=append_internal )
+      ENDIF
+   ENDIF
+
+END FUNCTION dom_def_att_real64
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Add attribute to database.
+!>
+!> @todo Try to combine similar code parts and shorten routine.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION dom_def_att_save( filename, variable_name, attribute, append ) RESULT( return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  filename                 !< name of file
+   CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  variable_name  !< name of variable
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_def_att_save'  !< name of routine
+
+   INTEGER(iwp) ::  a             !< loop index
+   INTEGER(iwp) ::  d             !< loop index
+   INTEGER(iwp) ::  f             !< loop index
+   INTEGER(iwp) ::  natt          !< number of attributes
+   INTEGER(iwp) ::  return_value  !< return value
+
+   LOGICAL             ::  found   !< true if variable or dimension of name 'variable_name' found
+   LOGICAL, INTENT(IN) ::  append  !< if true, append value to existing value
+
+   TYPE(attribute_type), INTENT(IN) ::  attribute  !< new attribute
+
+   TYPE(attribute_type), DIMENSION(:), ALLOCATABLE ::  atts_tmp  !< temporary attribute list
+
+
+   return_value = 0
+   found = .FALSE.
+
+   DO  f = 1, nf
+
+      IF ( TRIM( filename ) == files(f)%name )  THEN
+
+         !-- Add attribute to file
+         IF ( .NOT. PRESENT( variable_name ) )  THEN
+
+            !-- Initialize first file attribute
+            IF ( .NOT. ALLOCATED( files(f)%attributes ) )  THEN
+               natt = 1
+               ALLOCATE( files(f)%attributes(natt) )
+            ELSE
+               natt = SIZE( files(f)%attributes )
+
+               !-- Check if attribute already exists
+               DO  a = 1, natt
+                  IF ( files(f)%attributes(a)%name == attribute%name )  THEN
+                     IF ( append )  THEN
+                        !-- Append existing string attribute
+                        files(f)%attributes(a)%value_char =             &
+                           TRIM( files(f)%attributes(a)%value_char ) // &
+                           TRIM( attribute%value_char )
+                     ELSE
+                        files(f)%attributes(a) = attribute
+                     ENDIF
+                     found = .TRUE.
+                     EXIT
+                  ENDIF
+               ENDDO
+
+               !-- Extend attribute list by 1
+               IF ( .NOT. found )  THEN
+                  ALLOCATE( atts_tmp(natt) )
+                  atts_tmp = files(f)%attributes
+                  DEALLOCATE( files(f)%attributes )
+                  natt = natt + 1
+                  ALLOCATE( files(f)%attributes(natt) )
+                  files(f)%attributes(:natt-1) = atts_tmp
+                  DEALLOCATE( atts_tmp )
+               ENDIF
+            ENDIF
+
+            !-- Save new attribute to the end of the attribute list
+            IF ( .NOT. found )  THEN
+               files(f)%attributes(natt) = attribute
+               found = .TRUE.
+            ENDIF
+
+            EXIT
+
+         ELSE
+
+            !-- Add attribute to dimension
+            IF ( ALLOCATED( files(f)%dimensions ) )  THEN
+
+               DO  d = 1, SIZE( files(f)%dimensions )
+
+                  IF ( files(f)%dimensions(d)%name == TRIM( variable_name ) )  THEN
+
+                     IF ( .NOT. ALLOCATED( files(f)%dimensions(d)%attributes ) )  THEN
+                        !-- Initialize first attribute
+                        natt = 1
+                        ALLOCATE( files(f)%dimensions(d)%attributes(natt) )
+                     ELSE
+                        natt = SIZE( files(f)%dimensions(d)%attributes )
+
+                        !-- Check if attribute already exists
+                        DO  a = 1, natt
+                           IF ( files(f)%dimensions(d)%attributes(a)%name == attribute%name )  THEN
+                              IF ( append )  THEN
+                                 !-- Append existing character attribute
+                                 files(f)%dimensions(d)%attributes(a)%value_char =             &
+                                    TRIM( files(f)%dimensions(d)%attributes(a)%value_char ) // &
+                                    TRIM( attribute%value_char )
+                              ELSE
+                                 !-- Update existing attribute
+                                 files(f)%dimensions(d)%attributes(a) = attribute
+                              ENDIF
+                              found = .TRUE.
+                              EXIT
+                           ENDIF
+                        ENDDO
+
+                        !-- Extend attribute list
+                        IF ( .NOT. found )  THEN
+                           ALLOCATE( atts_tmp(natt) )
+                           atts_tmp = files(f)%dimensions(d)%attributes
+                           DEALLOCATE( files(f)%dimensions(d)%attributes )
+                           natt = natt + 1
+                           ALLOCATE( files(f)%dimensions(d)%attributes(natt) )
+                           files(f)%dimensions(d)%attributes(:natt-1) = atts_tmp
+                           DEALLOCATE( atts_tmp )
+                        ENDIF
+                     ENDIF
+
+                     !-- Add new attribute to database
+                     IF ( .NOT. found )  THEN
+                        files(f)%dimensions(d)%attributes(natt) = attribute
+                        found = .TRUE.
+                     ENDIF
+
+                     EXIT
+
+                  ENDIF  ! dimension found
+
+               ENDDO  ! loop over dimensions
+
+            ENDIF  ! dimensions exist in file
+
+            !-- Add attribute to variable
+            IF ( .NOT. found  .AND.  ALLOCATED( files(f)%variables) )  THEN
+
+               DO  d = 1, SIZE( files(f)%variables )
+
+                  IF ( files(f)%variables(d)%name == TRIM( variable_name ) )  THEN
+
+                     IF ( .NOT. ALLOCATED( files(f)%variables(d)%attributes ) )  THEN
+                        !-- Initialize first attribute
+                        natt = 1
+                        ALLOCATE( files(f)%variables(d)%attributes(natt) )
+                     ELSE
+                        natt = SIZE( files(f)%variables(d)%attributes )
+
+                        !-- Check if attribute already exists
+                        DO  a = 1, natt
+                           IF ( files(f)%variables(d)%attributes(a)%name == attribute%name )  THEN
+                              IF ( append )  THEN
+                                 !-- Append existing character attribute
+                                 files(f)%variables(d)%attributes(a)%value_char =             &
+                                    TRIM( files(f)%variables(d)%attributes(a)%value_char ) // &
+                                    TRIM( attribute%value_char )
+                              ELSE
+                                 !-- Update existing attribute
+                                 files(f)%variables(d)%attributes(a) = attribute
+                              ENDIF
+                              found = .TRUE.
+                              EXIT
+                           ENDIF
+                        ENDDO
+
+                        !-- Extend attribute list
+                        IF ( .NOT. found )  THEN
+                           ALLOCATE( atts_tmp(natt) )
+                           atts_tmp = files(f)%variables(d)%attributes
+                           DEALLOCATE( files(f)%variables(d)%attributes )
+                           natt = natt + 1
+                           ALLOCATE( files(f)%variables(d)%attributes(natt) )
+                           files(f)%variables(d)%attributes(:natt-1) = atts_tmp
+                           DEALLOCATE( atts_tmp )
+                        ENDIF
+
+                     ENDIF
+
+                     !-- Add new attribute to database
+                     IF ( .NOT. found )  THEN
+                        files(f)%variables(d)%attributes(natt) = attribute
+                        found = .TRUE.
+                     ENDIF
+
+                     EXIT
+
+                  ENDIF  ! variable found
+
+               ENDDO  ! loop over variables
+
+            ENDIF  ! variables exist in file
+
+            IF ( .NOT. found )  THEN
+               return_value = 1
+               CALL internal_message( 'error',                                        &
+                       routine_name //                                                &
+                       ': requested dimension/variable "' // TRIM( variable_name ) // &
+                       '" for attribute "' // TRIM( attribute%name ) //               &
+                       '" does not exist in file "' // TRIM( filename ) // '"' )
+            ENDIF
+
+            EXIT
+
+         ENDIF  ! variable_name present
+
+      ENDIF  ! check filename
+
+   ENDDO  ! loop over files
+
+   IF ( .NOT. found  .AND.  return_value == 0 )  THEN
+      return_value = 1
+      CALL internal_message( 'error',                                         &
+                             routine_name //                                  &
+                             ': requested file "' // TRIM( variable_name ) // &
+                             '" for attribute "' // TRIM( attribute%name ) // &
+                             '" does not exist' )
+   ENDIF
+
+END FUNCTION dom_def_att_save
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Start with output: clear database from unused files/dimensions, initialize
+!> files and write dimension values to files.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION dom_start_output() RESULT( return_value )
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_start_output'  !< name of routine
+
+   INTEGER(iwp) ::  d             !< loop index
+   INTEGER(iwp) ::  f             !< loop index
+   INTEGER(iwp) ::  return_value  !< return value
+
+   INTEGER(KIND=1), DIMENSION(:), ALLOCATABLE, TARGET ::  values_int8          !< target array for dimension values
+   INTEGER(KIND=1), DIMENSION(:), POINTER, CONTIGUOUS ::  values_int8_pointer  !< pointer to target array
+
+   INTEGER(KIND=2), DIMENSION(:), ALLOCATABLE, TARGET ::  values_int16          !< target array for dimension values
+   INTEGER(KIND=2), DIMENSION(:), POINTER, CONTIGUOUS ::  values_int16_pointer  !< pointer to target array
+
+   INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE, TARGET ::  values_int32          !< target array for dimension values
+   INTEGER(KIND=4), DIMENSION(:), POINTER, CONTIGUOUS ::  values_int32_pointer  !< pointer to target array
+
+   INTEGER(iwp), DIMENSION(:), ALLOCATABLE, TARGET ::  values_intwp          !< target array for dimension values
+   INTEGER(iwp), DIMENSION(:), POINTER, CONTIGUOUS ::  values_intwp_pointer  !< pointer to target array
+
+   REAL(KIND=4), DIMENSION(:), ALLOCATABLE, TARGET ::  values_real32          !< target array for dimension values
+   REAL(KIND=4), DIMENSION(:), POINTER, CONTIGUOUS ::  values_real32_pointer  !< pointer to target array
+
+   REAL(KIND=8), DIMENSION(:), ALLOCATABLE, TARGET ::  values_real64          !< target array for dimension values
+   REAL(KIND=8), DIMENSION(:), POINTER, CONTIGUOUS ::  values_real64_pointer  !< pointer to target array
+
+   REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET ::  values_realwp          !< target array for dimension values
+   REAL(wp), DIMENSION(:), POINTER, CONTIGUOUS ::  values_realwp_pointer  !< pointer to target array
+
+
+   !-- Clear database from empty files and unused dimensions
+   return_value = cleanup_database()
+
+   IF ( return_value == 0 )  THEN
+      DO  f = 1, nf
+
+         !-- Open file
+         CALL open_output_file( files(f)%format, files(f)%name, files(f)%id, &
+                                files(f)%is_init, return_value=return_value )
+
+         !-- Initialize file header:
+         !-- define dimensions and variables and write attributes
+         IF ( return_value == 0 )  &
+            CALL dom_init_file_header( files(f), return_value=return_value )
+
+         !-- End file definition
+         IF ( return_value == 0 )  &
+            CALL dom_init_end( files(f)%format, files(f)%id, return_value=return_value )
+
+         !-- Write dimension values into file
+         IF ( return_value == 0 )  THEN
+
+            DO  d = 1, SIZE( files(f)%dimensions )
+               IF ( ALLOCATED( files(f)%dimensions(d)%values_int8 ) )  THEN
+                  ALLOCATE( values_int8(files(f)%dimensions(d)%bounds(1): &
+                                        files(f)%dimensions(d)%bounds(2)) )
+                  values_int8 = files(f)%dimensions(d)%values_int8
+                  values_int8_pointer => values_int8
+                  return_value = dom_write_var( files(f)%name, files(f)%dimensions(d)%name, &
+                                    bounds_start=(/ files(f)%dimensions(d)%bounds(1) /),    &
+                                    bounds_end  =(/ files(f)%dimensions(d)%bounds(2) /),    &
+                                    var_int8_1d=values_int8_pointer )
+                  DEALLOCATE( values_int8 )
+               ELSEIF ( ALLOCATED( files(f)%dimensions(d)%values_int16 ) )  THEN
+                  ALLOCATE( values_int16(files(f)%dimensions(d)%bounds(1): &
+                                         files(f)%dimensions(d)%bounds(2)) )
+                  values_int16 = files(f)%dimensions(d)%values_int16
+                  values_int16_pointer => values_int16
+                  return_value = dom_write_var( files(f)%name, files(f)%dimensions(d)%name, &
+                                    bounds_start=(/ files(f)%dimensions(d)%bounds(1) /),    &
+                                    bounds_end  =(/ files(f)%dimensions(d)%bounds(2) /),    &
+                                    var_int16_1d=values_int16_pointer )
+                  DEALLOCATE( values_int16 )
+               ELSEIF ( ALLOCATED( files(f)%dimensions(d)%values_int32 ) )  THEN
+                  ALLOCATE( values_int32(files(f)%dimensions(d)%bounds(1): &
+                                         files(f)%dimensions(d)%bounds(2)) )
+                  values_int32 = files(f)%dimensions(d)%values_int32
+                  values_int32_pointer => values_int32
+                  return_value = dom_write_var( files(f)%name, files(f)%dimensions(d)%name, &
+                                    bounds_start=(/ files(f)%dimensions(d)%bounds(1) /),    &
+                                    bounds_end  =(/ files(f)%dimensions(d)%bounds(2) /),    &
+                                    var_int32_1d=values_int32_pointer )
+                  DEALLOCATE( values_int32 )
+               ELSEIF ( ALLOCATED( files(f)%dimensions(d)%values_intwp ) )  THEN
+                  ALLOCATE( values_intwp(files(f)%dimensions(d)%bounds(1): &
+                                         files(f)%dimensions(d)%bounds(2)) )
+                  values_intwp = files(f)%dimensions(d)%values_intwp
+                  values_intwp_pointer => values_intwp
+                  return_value = dom_write_var( files(f)%name, files(f)%dimensions(d)%name, &
+                                    bounds_start=(/ files(f)%dimensions(d)%bounds(1) /),    &
+                                    bounds_end  =(/ files(f)%dimensions(d)%bounds(2) /),    &
+                                    var_intwp_1d=values_intwp_pointer )
+                  DEALLOCATE( values_intwp )
+               ELSEIF ( ALLOCATED( files(f)%dimensions(d)%values_real32 ) )  THEN
+                  ALLOCATE( values_real32(files(f)%dimensions(d)%bounds(1): &
+                                          files(f)%dimensions(d)%bounds(2)) )
+                  values_real32 = files(f)%dimensions(d)%values_real32
+                  values_real32_pointer => values_real32
+                  return_value = dom_write_var( files(f)%name, files(f)%dimensions(d)%name, &
+                                    bounds_start=(/ files(f)%dimensions(d)%bounds(1) /),    &
+                                    bounds_end  =(/ files(f)%dimensions(d)%bounds(2) /),    &
+                                    var_real32_1d=values_real32_pointer )
+                  DEALLOCATE( values_real32 )
+               ELSEIF ( ALLOCATED( files(f)%dimensions(d)%values_real64 ) )  THEN
+                  ALLOCATE( values_real64(files(f)%dimensions(d)%bounds(1): &
+                                          files(f)%dimensions(d)%bounds(2)) )
+                  values_real64 = files(f)%dimensions(d)%values_real64
+                  values_real64_pointer => values_real64
+                  return_value = dom_write_var( files(f)%name, files(f)%dimensions(d)%name, &
+                                    bounds_start=(/ files(f)%dimensions(d)%bounds(1) /),    &
+                                    bounds_end  =(/ files(f)%dimensions(d)%bounds(2) /),    &
+                                    var_real64_1d=values_real64_pointer )
+                  DEALLOCATE( values_real64 )
+               ELSEIF ( ALLOCATED( files(f)%dimensions(d)%values_realwp ) )  THEN
+                  ALLOCATE( values_realwp(files(f)%dimensions(d)%bounds(1): &
+                                          files(f)%dimensions(d)%bounds(2)) )
+                  values_realwp = files(f)%dimensions(d)%values_realwp
+                  values_realwp_pointer => values_realwp
+                  return_value = dom_write_var( files(f)%name, files(f)%dimensions(d)%name, &
+                                    bounds_start=(/ files(f)%dimensions(d)%bounds(1) /),    &
+                                    bounds_end  =(/ files(f)%dimensions(d)%bounds(2) /),    &
+                                    var_realwp_1d=values_realwp_pointer )
+                  DEALLOCATE( values_realwp )
+               ENDIF
+               IF ( return_value /= 0 )  EXIT
+            ENDDO
+
+         ENDIF
+
+         IF ( return_value /= 0 )  EXIT
+
+      ENDDO
+   ENDIF
+
+   CALL internal_message( 'debug', routine_name // ': finished' )
+
+END FUNCTION dom_start_output
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check database and delete any unused dimensions and empty files (i.e. files
+!> without variables).
+!--------------------------------------------------------------------------------------------------!
+FUNCTION cleanup_database() RESULT( return_value )
+
+   ! CHARACTER(LEN=*), PARAMETER ::  routine_name = 'cleanup_database'  !< name of routine
+
+   INTEGER(iwp) ::  d             !< loop index
+   INTEGER(iwp) ::  f             !< loop index
+   INTEGER(iwp) ::  i             !< loop index
+   INTEGER(iwp) ::  ndim          !< number of dimensions in a file
+   INTEGER(iwp) ::  ndim_used     !< number of used dimensions in a file
+   INTEGER(iwp) ::  nf_used       !< number of used files
+   INTEGER(iwp) ::  nvar          !< number of variables in a file
+   INTEGER(iwp) ::  return_value  !< return value
+
+   LOGICAL, DIMENSION(1:nf)           ::  file_is_used       !< true if file contains variables
+   LOGICAL, DIMENSION(:), ALLOCATABLE ::  dimension_is_used  !< true if dimension is used by any variable
+
+   TYPE(dimension_type), DIMENSION(:), ALLOCATABLE ::  used_dimensions  !< list of used dimensions
+
+   TYPE(file_type), DIMENSION(:), ALLOCATABLE ::  used_files  !< list of used files
+
+
+   return_value = 0
+
+   !-- Set flag for files to be initialized
+   file_is_used(:) = .FALSE.
+   DO  f = 1, nf
+      IF ( ALLOCATED( files(f)%variables ) )  THEN
+         file_is_used(f) = .TRUE.
+      ENDIF
+   ENDDO
+
+   !-- Copy flagged files into temporary list
+   nf_used = COUNT( file_is_used )
+   ALLOCATE( used_files(nf_used) )
+   i = 0
+   DO  f = 1, nf
+      IF ( file_is_used(f) )  THEN
+         i = i + 1
+         used_files(i) = files(f)
+      ENDIF
+   ENDDO
+
+   !-- Replace file list with list of used files
+   DEALLOCATE( files )
+   nf = nf_used
+   ALLOCATE( files(nf) )
+   files = used_files
+   DEALLOCATE( used_files )
+
+   !-- Check every file for unused dimensions
+   DO  f = 1, nf
+
+      !-- Get number of defined dimensions
+      ndim = SIZE( files(f)%dimensions )
+      ALLOCATE( dimension_is_used(ndim) )
+
+      !-- Go through all variables and flag all used dimensions
+      nvar = SIZE( files(f)%variables )
+      DO  d = 1, ndim
+         DO  i = 1, nvar
+            dimension_is_used(d) = &
+               ANY( files(f)%dimensions(d)%name == files(f)%variables(i)%dimension_names )
+            IF ( dimension_is_used(d) )  EXIT
+         ENDDO
+      ENDDO
+
+      !-- Copy used dimensions to temporary list
+      ndim_used = COUNT( dimension_is_used )
+      ALLOCATE( used_dimensions(ndim_used) )
+      i = 0
+      DO  d = 1, ndim
+         IF ( dimension_is_used(d) )  THEN
+            i = i + 1
+            used_dimensions(i) = files(f)%dimensions(d)
+         ENDIF
+      ENDDO
+
+      !-- Replace dimension list with list of used dimensions
+      DEALLOCATE( files(f)%dimensions )
+      ndim = ndim_used
+      ALLOCATE( files(f)%dimensions(ndim) )
+      files(f)%dimensions = used_dimensions
+      DEALLOCATE( used_dimensions )
+      DEALLOCATE( dimension_is_used )
+
+   ENDDO
+
+END FUNCTION cleanup_database
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Open requested output file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE open_output_file( file_format, filename, file_id, is_init, return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  file_format  !< file format chosen for file
+   CHARACTER(LEN=*), INTENT(IN) ::  filename     !< name of file to be checked
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'open_output_file'  !< name of routine
+
+   INTEGER(iwp), INTENT(OUT) ::  file_id       !< file ID
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+   LOGICAL, INTENT(OUT) ::  is_init  !< true if file is opened
+
+
+   SELECT CASE ( TRIM( file_format ) )
+
+      CASE ( 'binary' )
+         CALL binary_open_file( filename, file_id, return_value )
+
+      CASE ( 'netcdf4-serial' )
+         CALL netcdf4_serial_open_file( filename, file_id, return_value )
+
+      CASE ( 'netcdf4-parallel' )
+         CALL netcdf4_parallel_open_file( filename, file_id, return_value )
+
+      CASE DEFAULT
+         return_value = 1
+         CALL internal_message( 'error', routine_name //                              &
+                                         ': file "' // TRIM( filename ) //            &
+                                         '": file format "' // TRIM( file_format ) // &
+                                         '" not supported' )
+
+   END SELECT
+
+   is_init = return_value == 0
+
+END SUBROUTINE open_output_file
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Define attributes, dimensions and variables.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE dom_init_file_header( file, return_value )
+
+   ! CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_init_file_header'  !< name of routine
+
+   INTEGER(iwp)              ::  a             !< loop index
+   INTEGER(iwp)              ::  d             !< loop index
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+   TYPE(file_type), INTENT(INOUT) ::  file  !< initialize header of this file
+
+
+   return_value  = 0
+
+   !-- Write file attributes
+   IF ( ALLOCATED( file%attributes ) )  THEN
+      DO  a = 1, SIZE( file%attributes )
+         return_value = write_attribute( file%format, file%id, var_id=no_var_id, &
+                                attribute=file%attributes(a) )
+         IF ( return_value /= 0 )  EXIT
+      ENDDO
+   ENDIF
+
+   IF ( return_value == 0 )  THEN
+
+      !-- Initialize file dimensions
+      DO  d = 1, SIZE( file%dimensions )
+
+         IF ( .NOT. file%dimensions(d)%is_masked )  THEN
+
+            !-- Initialize non-masked dimension
+            CALL init_file_dimension( file%format,                             &
+                    file%id, file%dimensions(d)%id, file%dimensions(d)%var_id, &
+                    file%dimensions(d)%name, file%dimensions(d)%data_type,     &
+                    file%dimensions(d)%length, file%dimensions(d)%is_init, return_value )
+
+         ELSE
+
+            !-- Initialize masked dimension
+            CALL init_file_dimension( file%format,                             &
+                    file%id, file%dimensions(d)%id, file%dimensions(d)%var_id, &
+                    file%dimensions(d)%name, file%dimensions(d)%data_type,     &
+                    file%dimensions(d)%length_mask, file%dimensions(d)%is_init, return_value )
+
+         ENDIF
+
+         IF ( return_value == 0  .AND.  ALLOCATED( file%dimensions(d)%attributes ) )  THEN
+            !-- Write dimension attributes
+            DO  a = 1, SIZE( file%dimensions(d)%attributes )
+               return_value = write_attribute( file%format, file%id, file%dimensions(d)%var_id, &
+                                 attribute=file%dimensions(d)%attributes(a) )
+               IF ( return_value /= 0 )  EXIT
+            ENDDO
+         ENDIF
+
+         IF ( return_value /= 0 )  EXIT
+
+      ENDDO
+
+      !-- Save dimension IDs for variables wihtin database
+      IF ( return_value == 0 )  &
+         CALL collect_dimesion_ids_for_variables( file%variables, file%dimensions, return_value )
+
+      !-- Initialize file variables
+      IF ( return_value == 0 )  THEN
+         DO  d = 1, SIZE( file%variables )
+
+            CALL init_file_variable( file%format, file%id, file%variables(d)%id, &
+                    file%variables(d)%name, file%variables(d)%data_type,         &
+                    file%variables(d)%dimension_ids,                             &
+                    file%variables(d)%is_init, file%variables(d)%is_global, return_value )
+
+            IF ( return_value == 0  .AND.  ALLOCATED( file%variables(d)%attributes ) )  THEN
+               !-- Write variable attribures
+               DO  a = 1, SIZE( file%variables(d)%attributes )
+                  return_value = write_attribute( file%format, file%id, file%variables(d)%id, &
+                                    attribute=file%variables(d)%attributes(a) )
+                  IF ( return_value /= 0 )  EXIT
+               ENDDO
+            ENDIF
+
+            IF ( return_value /= 0 )  EXIT
+
+         ENDDO
+      ENDIF
+
+   ENDIF
+
+END SUBROUTINE dom_init_file_header
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Write attribute to file.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION write_attribute( file_format, file_id, var_id, attribute ) RESULT( return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  file_format  !< file format chosen for file
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'write_attribute'  !< file format chosen for file
+
+   INTEGER(iwp) ::  file_id       !< file ID
+   INTEGER(iwp) ::  return_value  !< return value
+   INTEGER(iwp) ::  var_id        !< variable ID
+
+   TYPE(attribute_type), INTENT(IN) ::  attribute  !< attribute to be written
+
+
+   SELECT CASE ( TRIM( file_format ) )
+
+      CASE ( 'binary' )
+
+         SELECT CASE ( TRIM( attribute%data_type ) )
+
+            CASE( 'char' )
+               CALL binary_write_attribute( file_id=file_id, var_id=var_id,          &
+                       att_name=attribute%name, att_value_char=attribute%value_char, &
+                       return_value=return_value )
+
+            CASE( 'int8' )
+               CALL binary_write_attribute( file_id=file_id, var_id=var_id,          &
+                       att_name=attribute%name, att_value_int8=attribute%value_int8, &
+                       return_value=return_value )
+
+            CASE( 'int16' )
+               CALL binary_write_attribute( file_id=file_id, var_id=var_id,            &
+                       att_name=attribute%name, att_value_int16=attribute%value_int16, &
+                       return_value=return_value )
+
+            CASE( 'int32' )
+               CALL binary_write_attribute( file_id=file_id, var_id=var_id,            &
+                       att_name=attribute%name, att_value_int32=attribute%value_int32, &
+                       return_value=return_value )
+
+            CASE( 'real32' )
+               CALL binary_write_attribute( file_id=file_id, var_id=var_id,              &
+                       att_name=attribute%name, att_value_real32=attribute%value_real32, &
+                       return_value=return_value )
+
+            CASE( 'real64' )
+               CALL binary_write_attribute( file_id=file_id, var_id=var_id,              &
+                       att_name=attribute%name, att_value_real64=attribute%value_real64, &
+                       return_value=return_value )
+
+            CASE DEFAULT
+               return_value = 1
+               CALL internal_message( 'error',                                          &
+                                      routine_name //                                   &
+                                      ': attribute "' // TRIM( attribute%name ) //      &
+                                      '": data type "'// TRIM( attribute%data_type ) // &
+                                      '" not supported for file format "binary".' )
+
+         END SELECT
+
+      CASE ( 'netcdf4-serial' )
+
+         SELECT CASE ( TRIM( attribute%data_type ) )
+
+            CASE( 'char' )
+               CALL netcdf4_serial_write_attribute( file_id=file_id, var_id=var_id,  &
+                       att_name=attribute%name, att_value_char=attribute%value_char, &
+                       return_value=return_value )
+
+            CASE( 'int8' )
+               CALL netcdf4_serial_write_attribute( file_id=file_id, var_id=var_id,  &
+                       att_name=attribute%name, att_value_int8=attribute%value_int8, &
+                       return_value=return_value )
+
+            CASE( 'int16' )
+               CALL netcdf4_serial_write_attribute( file_id=file_id, var_id=var_id,    &
+                       att_name=attribute%name, att_value_int16=attribute%value_int16, &
+                       return_value=return_value )
+
+            CASE( 'int32' )
+               CALL netcdf4_serial_write_attribute( file_id=file_id, var_id=var_id,    &
+                       att_name=attribute%name, att_value_int32=attribute%value_int32, &
+                       return_value=return_value )
+
+            CASE( 'real32' )
+               CALL netcdf4_serial_write_attribute( file_id=file_id, var_id=var_id,      &
+                       att_name=attribute%name, att_value_real32=attribute%value_real32, &
+                       return_value=return_value )
+
+            CASE( 'real64' )
+               CALL netcdf4_serial_write_attribute( file_id=file_id, var_id=var_id,      &
+                       att_name=attribute%name, att_value_real64=attribute%value_real64, &
+                       return_value=return_value )
+
+            CASE DEFAULT
+               return_value = 1
+               CALL internal_message( 'error',                                          &
+                                      routine_name //                                   &
+                                      ': attribute "' // TRIM( attribute%name ) //      &
+                                      '": data type "'// TRIM( attribute%data_type ) // &
+                                      '" not supported for file format "netcdf4-serial".' )
+
+         END SELECT
+
+      CASE ( 'netcdf4-parallel' )
+
+         SELECT CASE ( TRIM( attribute%data_type ) )
+
+            CASE( 'char' )
+               CALL netcdf4_parallel_write_attribute( file_id=file_id, var_id=var_id, &
+                       att_name=attribute%name, att_value_char=attribute%value_char,  &
+                       return_value=return_value )
+
+            CASE( 'int8' )
+               CALL netcdf4_parallel_write_attribute( file_id=file_id, var_id=var_id, &
+                       att_name=attribute%name, att_value_int8=attribute%value_int8,  &
+                       return_value=return_value )
+
+            CASE( 'int16' )
+               CALL netcdf4_parallel_write_attribute( file_id=file_id, var_id=var_id,  &
+                       att_name=attribute%name, att_value_int16=attribute%value_int16, &
+                       return_value=return_value )
+
+            CASE( 'int32' )
+               CALL netcdf4_parallel_write_attribute( file_id=file_id, var_id=var_id,  &
+                       att_name=attribute%name, att_value_int32=attribute%value_int32, &
+                       return_value=return_value )
+
+            CASE( 'real32' )
+               CALL netcdf4_parallel_write_attribute( file_id=file_id, var_id=var_id,    &
+                       att_name=attribute%name, att_value_real32=attribute%value_real32, &
+                       return_value=return_value )
+
+            CASE( 'real64' )
+               CALL netcdf4_parallel_write_attribute( file_id=file_id, var_id=var_id,    &
+                       att_name=attribute%name, att_value_real64=attribute%value_real64, &
+                       return_value=return_value )
+
+            CASE DEFAULT
+               return_value = 1
+               CALL internal_message( 'error',                                          &
+                                      routine_name //                                   &
+                                      ': attribute "' // TRIM( attribute%name ) //      &
+                                      '": data type "'// TRIM( attribute%data_type ) // &
+                                      '" not supported for file format "netcdf4-parallel".' )
+
+         END SELECT
+
+      CASE DEFAULT
+         return_value = 1
+         CALL internal_message( 'error',        &
+                                routine_name // &
+                                ': unsupported file format "' // TRIM( file_format ) // '"' )
+
+   END SELECT
+
+END FUNCTION write_attribute
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize dimension in file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE init_file_dimension( file_format, file_id, dim_id, var_id,        &
+                                dim_name, dim_type, dim_length, is_init, return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  dim_name     !< name of dimension
+   CHARACTER(LEN=*), INTENT(IN) ::  dim_type     !< data type of dimension
+   CHARACTER(LEN=*), INTENT(IN) ::  file_format  !< file format chosen for file
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'init_file_dimension'  !< file format chosen for file
+
+   INTEGER(iwp), INTENT(OUT) ::  dim_id        !< dimension ID
+   INTEGER(iwp), INTENT(IN)  ::  dim_length    !< length of dimension
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+   INTEGER(iwp), INTENT(OUT) ::  var_id        !< associated variable ID
+
+   LOGICAL, INTENT(OUT) ::  is_init  !< true if dimension is initialized
+
+
+   SELECT CASE ( TRIM( file_format ) )
+
+      CASE ( 'binary' )
+         CALL binary_init_dimension( file_id, dim_id, var_id, &
+                 dim_name, dim_type, dim_length, is_init, return_value )
+
+      CASE ( 'netcdf4-serial' )
+         CALL netcdf4_serial_init_dimension( file_id, dim_id, var_id, &
+                 dim_name, dim_type, dim_length, is_init, return_value )
+
+      CASE ( 'netcdf4-parallel' )
+         CALL netcdf4_parallel_init_dimension( file_id, dim_id, var_id, &
+                 dim_name, dim_type, dim_length, is_init, return_value )
+
+      CASE DEFAULT
+         return_value = 1
+         WRITE( temp_string, * ) file_id
+         CALL internal_message( 'error', routine_name //                              &
+                                         ': file id = ' // TRIM( temp_string ) //     &
+                                         '": file format "' // TRIM( file_format ) // &
+                                         '" not supported' )
+
+   END SELECT
+
+END SUBROUTINE init_file_dimension
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Get dimension IDs and save them to variables.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE collect_dimesion_ids_for_variables( variables, dimensions, return_value )
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'collect_dimesion_ids_for_variables'  !< file format chosen for file
+
+   INTEGER(iwp) ::  d             !< loop index
+   INTEGER(iwp) ::  i             !< loop index
+   INTEGER(iwp) ::  j             !< loop index
+   INTEGER(iwp) ::  ndim          !< number of dimensions
+   INTEGER(iwp) ::  nvar          !< number of variables
+   INTEGER(iwp) ::  return_value  !< return value
+
+   LOGICAL ::  found  !< true if dimension required by variable was found in dimension list
+
+   TYPE(dimension_type), DIMENSION(:), INTENT(IN) ::  dimensions  !< list of dimensions in file
+
+   TYPE(variable_type), DIMENSION(:), INTENT(INOUT) ::  variables  !< list of variables in file
+
+
+   return_value  = 0
+   ndim = SIZE( dimensions )
+   nvar = SIZE( variables )
+
+   DO  i = 1, nvar
+      DO  j = 1, SIZE( variables(i)%dimension_names )
+         found = .FALSE.
+         DO  d = 1, ndim
+            IF ( variables(i)%dimension_names(j) == dimensions(d)%name )  THEN
+               variables(i)%dimension_ids(j) = dimensions(d)%id
+               found = .TRUE.
+               EXIT
+            ENDIF
+         ENDDO
+         IF ( .NOT. found )  THEN
+            return_value = 1
+            CALL internal_message( 'error',                                                 &
+                    routine_name // ': variable "' // TRIM( variables(i)%name ) //          &
+                    '": required dimension "' // TRIM( variables(i)%dimension_names(j) ) // &
+                    '" is undefined' )
+            EXIT
+         ENDIF
+      ENDDO
+      IF ( .NOT. found )  EXIT
+   ENDDO
+
+END SUBROUTINE collect_dimesion_ids_for_variables
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize variable.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE init_file_variable( file_format, file_id, var_id,   &
+                               var_name, var_type, var_dim_id, &
+                               is_init, is_global, return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  file_format  !< file format chosen for file
+   CHARACTER(LEN=*), INTENT(IN) ::  var_name     !< name of variable
+   CHARACTER(LEN=*), INTENT(IN) ::  var_type     !< data type of variable
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'init_file_variable'  !< file format chosen for file
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+   INTEGER(iwp), INTENT(OUT) ::  var_id        !< variable ID
+
+   INTEGER(iwp), DIMENSION(:), INTENT(IN) ::  var_dim_id  !< list of dimension IDs used by variable
+
+   LOGICAL, INTENT(IN)  ::  is_global  !< true if variable is global
+   LOGICAL, INTENT(OUT) ::  is_init    !< true if variable is initialized
+
+
+   SELECT CASE ( TRIM( file_format ) )
+
+      CASE ( 'binary' )
+         CALL binary_init_variable( file_id, var_id, var_name, var_type, var_dim_id, &
+                                    is_init, is_global, return_value )
+
+      CASE ( 'netcdf4-serial' )
+         CALL netcdf4_serial_init_variable( file_id, var_id, var_name, var_type, var_dim_id, &
+                                            is_init, is_global, return_value )
+
+      CASE ( 'netcdf4-parallel' )
+         CALL netcdf4_parallel_init_variable( file_id, var_id, var_name, var_type, var_dim_id, &
+                                              is_init, is_global, return_value )
+
+      CASE DEFAULT
+         return_value = 1
+         is_init = .FALSE.
+         CALL internal_message( 'error', routine_name // &
+                                ': unsupported file format "' // TRIM( file_format ) )
+
+   END SELECT
+
+END SUBROUTINE init_file_variable
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Finalize file definition/initialization.
+!>
+!> @todo Do we need an MPI barrier at the end?
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE dom_init_end( file_format, file_id, return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  file_format  !< file format
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_init_end'  !< name of routine
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file id
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+
+   SELECT CASE ( TRIM( file_format ) )
+
+      CASE ( 'binary' )
+         CALL binary_init_end( file_id, return_value )
+
+      CASE ( 'netcdf4-serial' )
+         CALL netcdf4_serial_init_end( file_id, return_value )
+
+      CASE ( 'netcdf4-parallel' )
+         CALL netcdf4_parallel_init_end( file_id, return_value )
+
+      CASE DEFAULT
+         return_value = 1
+         WRITE( temp_string, * ) file_id
+         CALL internal_message( 'error', routine_name //                             &
+                                         ': file id = ' // TRIM( temp_string ) //    &
+                                         ': file format "' // TRIM( file_format ) // &
+                                         '" not supported' )
+
+   END SELECT
+
+   ! CALL MPI_Barrier( MPI_COMM_WORLD, return_value )
+
+END SUBROUTINE dom_init_end
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Write variable to file.
+!> Example call:
+!>   dom_write_var( file_format = 'binary', &
+!>                  filename = 'DATA_OUTPUT_3D', &
+!>                  name = 'u', &
+!>                  var_real64_3d = u, &
+!>                  bounds_start = (/nzb, nys, nxl/), &
+!>                  bounds_end = (/nzt, nyn, nxr/)  )
+!--------------------------------------------------------------------------------------------------!
+FUNCTION dom_write_var( filename, name, bounds_start, bounds_end,       &
+            var_int8_0d,   var_int8_1d,   var_int8_2d,   var_int8_3d,   &
+            var_int16_0d,  var_int16_1d,  var_int16_2d,  var_int16_3d,  &
+            var_int32_0d,  var_int32_1d,  var_int32_2d,  var_int32_3d,  &
+            var_intwp_0d,  var_intwp_1d,  var_intwp_2d,  var_intwp_3d,  &
+            var_real32_0d, var_real32_1d, var_real32_2d, var_real32_3d, &
+            var_real64_0d, var_real64_1d, var_real64_2d, var_real64_3d, &
+            var_realwp_0d, var_realwp_1d, var_realwp_2d, var_realwp_3d  &
+            ) RESULT( return_value )
+
+   CHARACTER(LEN=charlen)       ::  file_format  !< file format chosen for file
+   CHARACTER(LEN=*), INTENT(IN) ::  filename     !< name of file
+   CHARACTER(LEN=*), INTENT(IN) ::  name         !< name of variable
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_write_var'  !< name of routine
+
+   INTEGER(iwp) ::  d             !< loop index
+   INTEGER(iwp) ::  file_id       !< file ID
+   INTEGER(iwp) ::  i             !< loop index
+   INTEGER(iwp) ::  j             !< loop index
+   INTEGER(iwp) ::  k             !< loop index
+   INTEGER(iwp) ::  return_value  !< return value
+   INTEGER(iwp) ::  var_id        !< variable ID
+
+   INTEGER(iwp), DIMENSION(:),   INTENT(IN)  ::  bounds_end        !< end index (upper bound) of variable at each dimension
+   INTEGER(iwp), DIMENSION(:),   INTENT(IN)  ::  bounds_start      !< start index (lower bound) of variable at each dimension
+   INTEGER(iwp), DIMENSION(:),   ALLOCATABLE ::  bounds_dim_start  !< start index (lower bound) of each dimension of variable
+   INTEGER(iwp), DIMENSION(:),   ALLOCATABLE ::  bounds_end_new    !< start index (upper bound) of masked variable at each dimension
+   INTEGER(iwp), DIMENSION(:),   ALLOCATABLE ::  bounds_start_new  !< start index (lower bound) of masked variable at each dimension
+   INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  masked_indices    !< dummy list holding all masked indices along a dimension
+
+   LOGICAL ::  do_output  !< true if any data lies within given range of masked dimension
+   LOGICAL ::  is_global  !< true if variable is global
+
+   INTEGER(KIND=1), POINTER, INTENT(IN), OPTIONAL                   ::  var_int8_0d  !< output variable
+   INTEGER(KIND=1), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_int8_1d  !< output variable
+   INTEGER(KIND=1), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_int8_2d  !< output variable
+   INTEGER(KIND=1), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_int8_3d  !< output variable
+
+   INTEGER(KIND=1), TARGET, ALLOCATABLE, DIMENSION(:)     ::  var_int8_1d_resorted  !< resorted output variable
+   INTEGER(KIND=1), TARGET, ALLOCATABLE, DIMENSION(:,:)   ::  var_int8_2d_resorted  !< resorted output variable
+   INTEGER(KIND=1), TARGET, ALLOCATABLE, DIMENSION(:,:,:) ::  var_int8_3d_resorted  !< resorted output variable
+
+   INTEGER(KIND=1), POINTER                               ::  var_int8_0d_pointer  !< output variable
+   INTEGER(KIND=1), POINTER, CONTIGUOUS, DIMENSION(:)     ::  var_int8_1d_pointer  !< output variable
+   INTEGER(KIND=1), POINTER, CONTIGUOUS, DIMENSION(:,:)   ::  var_int8_2d_pointer  !< output variable
+   INTEGER(KIND=1), POINTER, CONTIGUOUS, DIMENSION(:,:,:) ::  var_int8_3d_pointer  !< output variable
+
+   INTEGER(KIND=2), POINTER, INTENT(IN), OPTIONAL                   ::  var_int16_0d  !< output variable
+   INTEGER(KIND=2), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_int16_1d  !< output variable
+   INTEGER(KIND=2), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_int16_2d  !< output variable
+   INTEGER(KIND=2), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_int16_3d  !< output variable
+
+   INTEGER(KIND=2), TARGET, ALLOCATABLE, DIMENSION(:)     ::  var_int16_1d_resorted  !< resorted output variable
+   INTEGER(KIND=2), TARGET, ALLOCATABLE, DIMENSION(:,:)   ::  var_int16_2d_resorted  !< resorted output variable
+   INTEGER(KIND=2), TARGET, ALLOCATABLE, DIMENSION(:,:,:) ::  var_int16_3d_resorted  !< resorted output variable
+
+   INTEGER(KIND=2), POINTER                               ::  var_int16_0d_pointer  !< output variable
+   INTEGER(KIND=2), POINTER, CONTIGUOUS, DIMENSION(:)     ::  var_int16_1d_pointer  !< output variable
+   INTEGER(KIND=2), POINTER, CONTIGUOUS, DIMENSION(:,:)   ::  var_int16_2d_pointer  !< output variable
+   INTEGER(KIND=2), POINTER, CONTIGUOUS, DIMENSION(:,:,:) ::  var_int16_3d_pointer  !< output variable
+
+   INTEGER(KIND=4), POINTER, INTENT(IN), OPTIONAL                   ::  var_int32_0d  !< output variable
+   INTEGER(KIND=4), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_int32_1d  !< output variable
+   INTEGER(KIND=4), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_int32_2d  !< output variable
+   INTEGER(KIND=4), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_int32_3d  !< output variable
+
+   INTEGER(KIND=4), TARGET, ALLOCATABLE, DIMENSION(:)     ::  var_int32_1d_resorted  !< resorted output variable
+   INTEGER(KIND=4), TARGET, ALLOCATABLE, DIMENSION(:,:)   ::  var_int32_2d_resorted  !< resorted output variable
+   INTEGER(KIND=4), TARGET, ALLOCATABLE, DIMENSION(:,:,:) ::  var_int32_3d_resorted  !< resorted output variable
+
+   INTEGER(KIND=4), POINTER                               ::  var_int32_0d_pointer  !< output variable
+   INTEGER(KIND=4), POINTER, CONTIGUOUS, DIMENSION(:)     ::  var_int32_1d_pointer  !< output variable
+   INTEGER(KIND=4), POINTER, CONTIGUOUS, DIMENSION(:,:)   ::  var_int32_2d_pointer  !< output variable
+   INTEGER(KIND=4), POINTER, CONTIGUOUS, DIMENSION(:,:,:) ::  var_int32_3d_pointer  !< output variable
+
+   INTEGER(iwp), POINTER, INTENT(IN), OPTIONAL                   ::  var_intwp_0d  !< output variable
+   INTEGER(iwp), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_intwp_1d  !< output variable
+   INTEGER(iwp), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_intwp_2d  !< output variable
+   INTEGER(iwp), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_intwp_3d  !< output variable
+
+   INTEGER(iwp), TARGET, ALLOCATABLE, DIMENSION(:)     ::  var_intwp_1d_resorted  !< resorted output variable
+   INTEGER(iwp), TARGET, ALLOCATABLE, DIMENSION(:,:)   ::  var_intwp_2d_resorted  !< resorted output variable
+   INTEGER(iwp), TARGET, ALLOCATABLE, DIMENSION(:,:,:) ::  var_intwp_3d_resorted  !< resorted output variable
+
+   INTEGER(iwp), POINTER                               ::  var_intwp_0d_pointer  !< output variable
+   INTEGER(iwp), POINTER, CONTIGUOUS, DIMENSION(:)     ::  var_intwp_1d_pointer  !< output variable
+   INTEGER(iwp), POINTER, CONTIGUOUS, DIMENSION(:,:)   ::  var_intwp_2d_pointer  !< output variable
+   INTEGER(iwp), POINTER, CONTIGUOUS, DIMENSION(:,:,:) ::  var_intwp_3d_pointer  !< output variable
+
+   REAL(KIND=4), POINTER, INTENT(IN), OPTIONAL                   ::  var_real32_0d  !< output variable
+   REAL(KIND=4), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_real32_1d  !< output variable
+   REAL(KIND=4), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_real32_2d  !< output variable
+   REAL(KIND=4), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_real32_3d  !< output variable
+
+   REAL(KIND=4), TARGET, ALLOCATABLE, DIMENSION(:)     ::  var_real32_1d_resorted  !< resorted output variable
+   REAL(KIND=4), TARGET, ALLOCATABLE, DIMENSION(:,:)   ::  var_real32_2d_resorted  !< resorted output variable
+   REAL(KIND=4), TARGET, ALLOCATABLE, DIMENSION(:,:,:) ::  var_real32_3d_resorted  !< resorted output variable
+
+   REAL(KIND=4), POINTER                               ::  var_real32_0d_pointer  !< output variable
+   REAL(KIND=4), POINTER, CONTIGUOUS, DIMENSION(:)     ::  var_real32_1d_pointer  !< output variable
+   REAL(KIND=4), POINTER, CONTIGUOUS, DIMENSION(:,:)   ::  var_real32_2d_pointer  !< output variable
+   REAL(KIND=4), POINTER, CONTIGUOUS, DIMENSION(:,:,:) ::  var_real32_3d_pointer  !< output variable
+
+   REAL(KIND=8), POINTER, INTENT(IN), OPTIONAL                   ::  var_real64_0d  !< output variable
+   REAL(KIND=8), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_real64_1d  !< output variable
+   REAL(KIND=8), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_real64_2d  !< output variable
+   REAL(KIND=8), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_real64_3d  !< output variable
+
+   REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:)     ::  var_real64_1d_resorted  !< resorted output variable
+   REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:,:)   ::  var_real64_2d_resorted  !< resorted output variable
+   REAL(KIND=8), TARGET, ALLOCATABLE, DIMENSION(:,:,:) ::  var_real64_3d_resorted  !< resorted output variable
+
+   REAL(KIND=8), POINTER                               ::  var_real64_0d_pointer  !< output variable
+   REAL(KIND=8), POINTER, CONTIGUOUS, DIMENSION(:)     ::  var_real64_1d_pointer  !< output variable
+   REAL(KIND=8), POINTER, CONTIGUOUS, DIMENSION(:,:)   ::  var_real64_2d_pointer  !< output variable
+   REAL(KIND=8), POINTER, CONTIGUOUS, DIMENSION(:,:,:) ::  var_real64_3d_pointer  !< output variable
+
+   REAL(wp), POINTER, INTENT(IN), OPTIONAL                   ::  var_realwp_0d  !< output variable
+   REAL(wp), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_realwp_1d  !< output variable
+   REAL(wp), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_realwp_2d  !< output variable
+   REAL(wp), POINTER, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_realwp_3d  !< output variable
+
+   REAL(wp), TARGET, ALLOCATABLE, DIMENSION(:)     ::  var_realwp_1d_resorted  !< resorted output variable
+   REAL(wp), TARGET, ALLOCATABLE, DIMENSION(:,:)   ::  var_realwp_2d_resorted  !< resorted output variable
+   REAL(wp), TARGET, ALLOCATABLE, DIMENSION(:,:,:) ::  var_realwp_3d_resorted  !< resorted output variable
+
+   REAL(wp), POINTER                               ::  var_realwp_0d_pointer  !< output variable
+   REAL(wp), POINTER, CONTIGUOUS, DIMENSION(:)     ::  var_realwp_1d_pointer  !< output variable
+   REAL(wp), POINTER, CONTIGUOUS, DIMENSION(:,:)   ::  var_realwp_2d_pointer  !< output variable
+   REAL(wp), POINTER, CONTIGUOUS, DIMENSION(:,:,:) ::  var_realwp_3d_pointer  !< output variable
+
+   TYPE(dimension_type), DIMENSION(:), ALLOCATABLE ::  dimension_list  !< list of used dimensions of variable
+
+
+   CALL internal_message( 'debug', routine_name // ': write ' // TRIM( name ) // &
+                                   ' into file ' // TRIM( filename ) )
+
+   !-- Search for variable within file
+   CALL find_var_in_file( filename, name, file_format, file_id, var_id, &
+                      is_global, dimension_list, return_value=return_value  )
+
+   IF ( return_value == 0 )  THEN
+
+      !-- Check if the correct amount of variable bounds were given
+      IF ( SIZE( bounds_start ) /= SIZE( dimension_list )  .OR.  &
+           SIZE( bounds_end ) /= SIZE( dimension_list ) )  THEN
+         return_value = 1
+         CALL internal_message( 'error', routine_name //                      &
+                                         ': variable "' // TRIM( name ) //    &
+                                         '" in file "' // TRIM( filename ) // &
+                                         '": given bounds do not match with number of dimensions' )
+      ENDIF
+
+   ENDIF
+
+
+   IF ( return_value == 0 )  THEN
+
+      !-- Save starting index (lower bounds) of each dimension
+      ALLOCATE( bounds_dim_start(SIZE( dimension_list )) )
+      DO  d = 1, SIZE( dimension_list )
+         bounds_dim_start(d) = dimension_list(d)%bounds(1)
+      ENDDO
+
+      ALLOCATE( bounds_start_new(SIZE( dimension_list )) )
+      ALLOCATE( bounds_end_new(SIZE( dimension_list ))   )
+
+      WRITE( temp_string, * ) bounds_start
+      CALL internal_message( 'debug', routine_name //                     &
+                                      ': file "' // TRIM( filename ) //   &
+                                       '": variable "' // TRIM( name ) // &
+                                       '": bounds_start =' // TRIM( temp_string ) )
+      WRITE( temp_string, * ) bounds_end
+      CALL internal_message( 'debug', routine_name //                     &
+                                      ': file "' // TRIM( filename ) //   &
+                                       '": variable "' // TRIM( name ) // &
+                                       '": bounds_end =' // TRIM( temp_string ) )
+
+      !-- Get bounds for masking
+      CALL get_masked_indices_and_masked_dimension_bounds( dimension_list, &
+              bounds_start, bounds_end, bounds_start_new, bounds_end_new,  &
+              masked_indices, do_output )
+
+      WRITE( temp_string, * ) bounds_start_new
+      CALL internal_message( 'debug', routine_name //                     &
+                                      ': file "' // TRIM( filename ) //   &
+                                       '": variable "' // TRIM( name ) // &
+                                       '": bounds_start_new =' // TRIM( temp_string ) )
+      WRITE( temp_string, * ) bounds_end_new
+      CALL internal_message( 'debug', routine_name //                     &
+                                      ': file "' // TRIM( filename ) //   &
+                                       '": variable "' // TRIM( name ) // &
+                                       '": bounds_end_new =' // TRIM( temp_string ) )
+
+      !-- Mask and resort variable
+      !-- 8bit integer output
+      IF ( PRESENT( var_int8_0d ) )  THEN
+         var_int8_0d_pointer => var_int8_0d
+      ELSEIF ( PRESENT( var_int8_1d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_int8_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            !$OMP PARALLEL PRIVATE (i)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               var_int8_1d_resorted(i) = var_int8_1d(masked_indices(1,i))
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_int8_1d_resorted(1) )
+            var_int8_1d_resorted = 0_1
+         ENDIF
+         var_int8_1d_pointer => var_int8_1d_resorted
+      ELSEIF ( PRESENT( var_int8_2d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_int8_2d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                           bounds_start_new(2):bounds_end_new(2)) )
+            !$OMP PARALLEL PRIVATE (i,j)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  var_int8_2d_resorted(i,j) = var_int8_2d(masked_indices(2,j), &
+                                                          masked_indices(1,i)  )
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_int8_2d_resorted(1,1) )
+            var_int8_2d_resorted = 0_1
+         ENDIF
+         var_int8_2d_pointer => var_int8_2d_resorted
+      ELSEIF ( PRESENT( var_int8_3d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_int8_3d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                           bounds_start_new(2):bounds_end_new(2), &
+                                           bounds_start_new(3):bounds_end_new(3)) )
+            !$OMP PARALLEL PRIVATE (i,j,k)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  DO  k = bounds_start_new(3), bounds_end_new(3)
+                     var_int8_3d_resorted(i,j,k) = var_int8_3d(masked_indices(3,k), &
+                                                               masked_indices(2,j), &
+                                                               masked_indices(1,i)  )
+                  ENDDO
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_int8_3d_resorted(1,1,1) )
+            var_int8_3d_resorted = 0_1
+         ENDIF
+         var_int8_3d_pointer => var_int8_3d_resorted
+
+      !-- 16bit integer output
+      ELSEIF ( PRESENT( var_int16_0d ) )  THEN
+         var_int16_0d_pointer => var_int16_0d
+      ELSEIF ( PRESENT( var_int16_1d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_int16_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            !$OMP PARALLEL PRIVATE (i)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               var_int16_1d_resorted(i) = var_int16_1d(masked_indices(1,i))
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_int16_1d_resorted(1) )
+            var_int16_1d_resorted = 0_1
+         ENDIF
+         var_int16_1d_pointer => var_int16_1d_resorted
+      ELSEIF ( PRESENT( var_int16_2d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_int16_2d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                            bounds_start_new(2):bounds_end_new(2)) )
+            !$OMP PARALLEL PRIVATE (i,j)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  var_int16_2d_resorted(i,j) = var_int16_2d(masked_indices(2,j), &
+                                                            masked_indices(1,i))
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_int16_2d_resorted(1,1) )
+            var_int16_2d_resorted = 0_1
+         ENDIF
+         var_int16_2d_pointer => var_int16_2d_resorted
+      ELSEIF ( PRESENT( var_int16_3d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_int16_3d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                            bounds_start_new(2):bounds_end_new(2), &
+                                            bounds_start_new(3):bounds_end_new(3)) )
+            !$OMP PARALLEL PRIVATE (i,j,k)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  DO  k = bounds_start_new(3), bounds_end_new(3)
+                     var_int16_3d_resorted(i,j,k) = var_int16_3d(masked_indices(3,k), &
+                                                                 masked_indices(2,j), &
+                                                                 masked_indices(1,i)  )
+                  ENDDO
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_int16_3d_resorted(1,1,1) )
+            var_int16_3d_resorted = 0_1
+         ENDIF
+         var_int16_3d_pointer => var_int16_3d_resorted
+
+      !-- 32bit integer output
+      ELSEIF ( PRESENT( var_int32_0d ) )  THEN
+         var_int32_0d_pointer => var_int32_0d
+      ELSEIF ( PRESENT( var_int32_1d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_int32_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            !$OMP PARALLEL PRIVATE (i)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               var_int32_1d_resorted(i) = var_int32_1d(masked_indices(1,i))
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_int32_1d_resorted(1) )
+            var_int32_1d_resorted = 0_1
+         ENDIF
+         var_int32_1d_pointer => var_int32_1d_resorted
+      ELSEIF ( PRESENT( var_int32_2d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_int32_2d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                            bounds_start_new(2):bounds_end_new(2)) )
+            !$OMP PARALLEL PRIVATE (i,j)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  var_int32_2d_resorted(i,j) = var_int32_2d(masked_indices(2,j), &
+                                                            masked_indices(1,i)  )
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_int32_2d_resorted(1,1) )
+            var_int32_2d_resorted = 0_1
+         ENDIF
+         var_int32_2d_pointer => var_int32_2d_resorted
+      ELSEIF ( PRESENT( var_int32_3d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_int32_3d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                            bounds_start_new(2):bounds_end_new(2), &
+                                            bounds_start_new(3):bounds_end_new(3)) )
+            !$OMP PARALLEL PRIVATE (i,j,k)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  DO  k = bounds_start_new(3), bounds_end_new(3)
+                     var_int32_3d_resorted(i,j,k) = var_int32_3d(masked_indices(3,k), &
+                                                                 masked_indices(2,j), &
+                                                                 masked_indices(1,i)  )
+                  ENDDO
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_int32_3d_resorted(1,1,1) )
+            var_int32_3d_resorted = 0_1
+         ENDIF
+         var_int32_3d_pointer => var_int32_3d_resorted
+
+      !-- working-precision integer output
+      ELSEIF ( PRESENT( var_intwp_0d ) )  THEN
+         var_intwp_0d_pointer => var_intwp_0d
+      ELSEIF ( PRESENT( var_intwp_1d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_intwp_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            !$OMP PARALLEL PRIVATE (i)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               var_intwp_1d_resorted(i) = var_intwp_1d(masked_indices(1,i))
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_intwp_1d_resorted(1) )
+            var_intwp_1d_resorted = 0_1
+         ENDIF
+         var_intwp_1d_pointer => var_intwp_1d_resorted
+      ELSEIF ( PRESENT( var_intwp_2d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_intwp_2d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                            bounds_start_new(2):bounds_end_new(2)) )
+            !$OMP PARALLEL PRIVATE (i,j)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  var_intwp_2d_resorted(i,j) = var_intwp_2d(masked_indices(2,j), &
+                                                            masked_indices(1,i)  )
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_intwp_2d_resorted(1,1) )
+            var_intwp_2d_resorted = 0_1
+         ENDIF
+         var_intwp_2d_pointer => var_intwp_2d_resorted
+      ELSEIF ( PRESENT( var_intwp_3d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_intwp_3d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                            bounds_start_new(2):bounds_end_new(2), &
+                                            bounds_start_new(3):bounds_end_new(3)) )
+            !$OMP PARALLEL PRIVATE (i,j,k)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  DO  k = bounds_start_new(3), bounds_end_new(3)
+                     var_intwp_3d_resorted(i,j,k) = var_intwp_3d(masked_indices(3,k), &
+                                                                 masked_indices(2,j), &
+                                                                 masked_indices(1,i)  )
+                  ENDDO
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_intwp_3d_resorted(1,1,1) )
+            var_intwp_3d_resorted = 0_1
+         ENDIF
+         var_intwp_3d_pointer => var_intwp_3d_resorted
+
+      !-- 32bit real output
+      ELSEIF ( PRESENT( var_real32_0d ) )  THEN
+         var_real32_0d_pointer => var_real32_0d
+      ELSEIF ( PRESENT( var_real32_1d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_real32_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            !$OMP PARALLEL PRIVATE (i)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               var_real32_1d_resorted(i) = var_real32_1d(masked_indices(1,i))
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_real32_1d_resorted(1) )
+            var_real32_1d_resorted = 0_1
+         ENDIF
+         var_real32_1d_pointer => var_real32_1d_resorted
+      ELSEIF ( PRESENT( var_real32_2d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_real32_2d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                             bounds_start_new(2):bounds_end_new(2)) )
+            !$OMP PARALLEL PRIVATE (i,j)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  var_real32_2d_resorted(i,j) = var_real32_2d(masked_indices(2,j), &
+                                                              masked_indices(1,i)  )
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_real32_2d_resorted(1,1) )
+            var_real32_2d_resorted = 0_1
+         ENDIF
+         var_real32_2d_pointer => var_real32_2d_resorted
+      ELSEIF ( PRESENT( var_real32_3d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_real32_3d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                             bounds_start_new(2):bounds_end_new(2), &
+                                             bounds_start_new(3):bounds_end_new(3)) )
+            !$OMP PARALLEL PRIVATE (i,j,k)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  DO  k = bounds_start_new(3), bounds_end_new(3)
+                     var_real32_3d_resorted(i,j,k) = var_real32_3d(masked_indices(3,k), &
+                                                                   masked_indices(2,j), &
+                                                                   masked_indices(1,i)  )
+                  ENDDO
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_real32_3d_resorted(1,1,1) )
+            var_real32_3d_resorted = 0_1
+         ENDIF
+         var_real32_3d_pointer => var_real32_3d_resorted
+
+      !-- 64bit real output
+      ELSEIF ( PRESENT( var_real64_0d ) )  THEN
+         var_real64_0d_pointer => var_real64_0d
+      ELSEIF ( PRESENT( var_real64_1d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_real64_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            !$OMP PARALLEL PRIVATE (i)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               var_real64_1d_resorted(i) = var_real64_1d(masked_indices(1,i))
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_real64_1d_resorted(1) )
+            var_real64_1d_resorted = 0_1
+         ENDIF
+         var_real64_1d_pointer => var_real64_1d_resorted
+      ELSEIF ( PRESENT( var_real64_2d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_real64_2d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                             bounds_start_new(2):bounds_end_new(2)) )
+            !$OMP PARALLEL PRIVATE (i,j)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  var_real64_2d_resorted(i,j) = var_real64_2d(masked_indices(2,j), &
+                                                              masked_indices(1,i)  )
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_real64_2d_resorted(1,1) )
+            var_real64_2d_resorted = 0_1
+         ENDIF
+         var_real64_2d_pointer => var_real64_2d_resorted
+      ELSEIF ( PRESENT( var_real64_3d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_real64_3d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                             bounds_start_new(2):bounds_end_new(2), &
+                                             bounds_start_new(3):bounds_end_new(3)) )
+            !$OMP PARALLEL PRIVATE (i,j,k)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  DO  k = bounds_start_new(3), bounds_end_new(3)
+                     var_real64_3d_resorted(i,j,k) = var_real64_3d(masked_indices(3,k), &
+                                                                   masked_indices(2,j), &
+                                                                   masked_indices(1,i)  )
+                  ENDDO
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_real64_3d_resorted(1,1,1) )
+            var_real64_3d_resorted = 0_1
+         ENDIF
+         var_real64_3d_pointer => var_real64_3d_resorted
+
+      !-- working-precision real output
+      ELSEIF ( PRESENT( var_realwp_0d ) )  THEN
+         var_realwp_0d_pointer => var_realwp_0d
+      ELSEIF ( PRESENT( var_realwp_1d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_realwp_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            !$OMP PARALLEL PRIVATE (i)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               var_realwp_1d_resorted(i) = var_realwp_1d(masked_indices(1,i))
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_realwp_1d_resorted(1) )
+            var_realwp_1d_resorted = 0_1
+         ENDIF
+         var_realwp_1d_pointer => var_realwp_1d_resorted
+      ELSEIF ( PRESENT( var_realwp_2d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_realwp_2d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                             bounds_start_new(2):bounds_end_new(2)) )
+            !$OMP PARALLEL PRIVATE (i,j)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  var_realwp_2d_resorted(i,j) = var_realwp_2d(masked_indices(2,j), &
+                                                              masked_indices(1,i)  )
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_realwp_2d_resorted(1,1) )
+            var_realwp_2d_resorted = 0_1
+         ENDIF
+         var_realwp_2d_pointer => var_realwp_2d_resorted
+      ELSEIF ( PRESENT( var_realwp_3d ) )  THEN
+         IF ( do_output ) THEN
+            ALLOCATE( var_realwp_3d_resorted(bounds_start_new(1):bounds_end_new(1), &
+                                             bounds_start_new(2):bounds_end_new(2), &
+                                             bounds_start_new(3):bounds_end_new(3)) )
+            !$OMP PARALLEL PRIVATE (i,j,k)
+            !$OMP DO
+            DO  i = bounds_start_new(1), bounds_end_new(1)
+               DO  j = bounds_start_new(2), bounds_end_new(2)
+                  DO  k = bounds_start_new(3), bounds_end_new(3)
+                     var_realwp_3d_resorted(i,j,k) = var_realwp_3d(masked_indices(3,k), &
+                                                                   masked_indices(2,j), &
+                                                                   masked_indices(1,i)  )
+                  ENDDO
+               ENDDO
+            ENDDO
+            !$OMP END PARALLEL
+         ELSE
+            ALLOCATE( var_realwp_3d_resorted(1,1,1) )
+            var_realwp_3d_resorted = 0_1
+         ENDIF
+         var_realwp_3d_pointer => var_realwp_3d_resorted
+
+      ELSE
+         return_value = 1
+         CALL internal_message( 'error', routine_name //                      &
+                                         ': variable "' // TRIM( name ) //    &
+                                         '" in file "' // TRIM( filename ) // &
+                                         '": no values given to output' )
+      ENDIF
+
+      DEALLOCATE( masked_indices )
+
+   ENDIF  ! Check for error
+
+   IF ( return_value == 0 )  THEN
+
+      !-- Write variable into file
+      SELECT CASE ( TRIM( file_format ) )
+
+         CASE ( 'binary' )
+            !-- 8bit integer output
+            IF ( PRESENT( var_int8_0d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int8_0d=var_int8_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int8_1d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int8_1d=var_int8_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int8_2d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int8_2d=var_int8_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int8_3d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int8_3d=var_int8_3d_pointer, return_value=return_value )
+            !-- 16bit integer output
+            ELSEIF ( PRESENT( var_int16_0d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int16_0d=var_int16_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int16_1d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int16_1d=var_int16_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int16_2d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int16_2d=var_int16_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int16_3d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int16_3d=var_int16_3d_pointer, return_value=return_value )
+            !-- 32bit integer output
+            ELSEIF ( PRESENT( var_int32_0d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int32_0d=var_int32_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int32_1d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int32_1d=var_int32_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int32_2d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int32_2d=var_int32_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int32_3d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int32_3d=var_int32_3d_pointer, return_value=return_value )
+            !-- working-precision integer output
+            ELSEIF ( PRESENT( var_intwp_0d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_intwp_0d=var_intwp_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_intwp_1d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_intwp_1d=var_intwp_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_intwp_2d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_intwp_2d=var_intwp_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_intwp_3d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_intwp_3d=var_intwp_3d_pointer, return_value=return_value )
+            !-- 32bit real output
+            ELSEIF ( PRESENT( var_real32_0d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real32_0d=var_real32_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real32_1d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real32_1d=var_real32_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real32_2d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real32_2d=var_real32_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real32_3d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real32_3d=var_real32_3d_pointer, return_value=return_value )
+            !-- 64bit real output
+            ELSEIF ( PRESENT( var_real64_0d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real64_0d=var_real64_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real64_1d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real64_1d=var_real64_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real64_2d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real64_2d=var_real64_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real64_3d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real64_3d=var_real64_3d_pointer, return_value=return_value )
+            !-- working-precision real output
+            ELSEIF ( PRESENT( var_realwp_0d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_realwp_0d=var_realwp_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_realwp_1d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_realwp_1d=var_realwp_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_realwp_2d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_realwp_2d=var_realwp_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_realwp_3d ) )  THEN
+               CALL binary_write_variable( file_id, var_id,                                      &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_realwp_3d=var_realwp_3d_pointer, return_value=return_value )
+            ELSE
+               return_value = 1
+               CALL internal_message( 'error', routine_name //                                    &
+                                               ': variable "' // TRIM( name ) //                  &
+                                               '" in file "' // TRIM( filename ) //               &
+                                               '": output_type not supported by file format "' // &
+                                               TRIM( file_format ) // '"' )
+            ENDIF
+
+         CASE ( 'netcdf4-serial' )
+            !-- 8bit integer output
+            IF ( PRESENT( var_int8_0d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int8_0d=var_int8_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int8_1d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int8_1d=var_int8_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int8_2d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int8_2d=var_int8_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int8_3d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int8_3d=var_int8_3d_pointer, return_value=return_value )
+            !-- 16bit integer output
+            ELSEIF ( PRESENT( var_int16_0d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int16_0d=var_int16_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int16_1d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int16_1d=var_int16_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int16_2d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int16_2d=var_int16_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int16_3d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int16_3d=var_int16_3d_pointer, return_value=return_value )
+            !-- 32bit integer output
+            ELSEIF ( PRESENT( var_int32_0d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int32_0d=var_int32_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int32_1d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int32_1d=var_int32_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int32_2d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int32_2d=var_int32_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int32_3d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int32_3d=var_int32_3d_pointer, return_value=return_value )
+            !-- working-precision integer output
+            ELSEIF ( PRESENT( var_intwp_0d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_intwp_0d=var_intwp_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_intwp_1d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_intwp_1d=var_intwp_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_intwp_2d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_intwp_2d=var_intwp_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_intwp_3d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_intwp_3d=var_intwp_3d_pointer, return_value=return_value )
+            !-- 32bit real output
+            ELSEIF ( PRESENT( var_real32_0d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real32_0d=var_real32_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real32_1d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real32_1d=var_real32_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real32_2d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real32_2d=var_real32_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real32_3d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real32_3d=var_real32_3d_pointer, return_value=return_value )
+            !-- 64bit real output
+            ELSEIF ( PRESENT( var_real64_0d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real64_0d=var_real64_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real64_1d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real64_1d=var_real64_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real64_2d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real64_2d=var_real64_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real64_3d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real64_3d=var_real64_3d_pointer, return_value=return_value )
+            !-- working-precision real output
+            ELSEIF ( PRESENT( var_realwp_0d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_realwp_0d=var_realwp_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_realwp_1d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_realwp_1d=var_realwp_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_realwp_2d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_realwp_2d=var_realwp_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_realwp_3d ) )  THEN
+               CALL netcdf4_serial_write_variable( file_id, var_id,                              &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_realwp_3d=var_realwp_3d_pointer, return_value=return_value )
+            ELSE
+               return_value = 1
+               CALL internal_message( 'error', routine_name //                                    &
+                                               ': variable "' // TRIM( name ) //                  &
+                                               '" in file "' // TRIM( filename ) //               &
+                                               '": output_type not supported by file format "' // &
+                                               TRIM( file_format ) // '"' )
+            ENDIF
+
+         CASE ( 'netcdf4-parallel' )
+            !-- 8bit integer output
+            IF ( PRESENT( var_int8_0d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int8_0d=var_int8_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int8_1d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int8_1d=var_int8_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int8_2d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int8_2d=var_int8_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int8_3d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int8_3d=var_int8_3d_pointer, return_value=return_value )
+            !-- 16bit integer output
+            ELSEIF ( PRESENT( var_int16_0d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int16_0d=var_int16_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int16_1d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int16_1d=var_int16_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int16_2d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int16_2d=var_int16_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int16_3d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int16_3d=var_int16_3d_pointer, return_value=return_value )
+            !-- 32bit integer output
+            ELSEIF ( PRESENT( var_int32_0d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int32_0d=var_int32_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int32_1d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int32_1d=var_int32_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int32_2d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int32_2d=var_int32_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_int32_3d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_int32_3d=var_int32_3d_pointer, return_value=return_value )
+            !-- working-precision integer output
+            ELSEIF ( PRESENT( var_intwp_0d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_intwp_0d=var_intwp_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_intwp_1d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_intwp_1d=var_intwp_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_intwp_2d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_intwp_2d=var_intwp_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_intwp_3d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_intwp_3d=var_intwp_3d_pointer, return_value=return_value )
+            !-- 32bit real output
+            ELSEIF ( PRESENT( var_real32_0d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real32_0d=var_real32_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real32_1d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real32_1d=var_real32_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real32_2d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real32_2d=var_real32_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real32_3d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real32_3d=var_real32_3d_pointer, return_value=return_value )
+            !-- 64bit real output
+            ELSEIF ( PRESENT( var_real64_0d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real64_0d=var_real64_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real64_1d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real64_1d=var_real64_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real64_2d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real64_2d=var_real64_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_real64_3d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_real64_3d=var_real64_3d_pointer, return_value=return_value )
+            !-- working-precision real output
+            ELSEIF ( PRESENT( var_realwp_0d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_realwp_0d=var_realwp_0d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_realwp_1d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_realwp_1d=var_realwp_1d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_realwp_2d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_realwp_2d=var_realwp_2d_pointer, return_value=return_value )
+            ELSEIF ( PRESENT( var_realwp_3d ) )  THEN
+               CALL netcdf4_parallel_write_variable( file_id, var_id,                            &
+                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+                       var_realwp_3d=var_realwp_3d_pointer, return_value=return_value )
+            ELSE
+               return_value = 1
+               CALL internal_message( 'error', routine_name //                                    &
+                                               ': variable "' // TRIM( name ) //                  &
+                                               '" in file "' // TRIM( filename ) //               &
+                                               '": output_type not supported by file format "' // &
+                                               TRIM( file_format ) // '"' )
+            ENDIF
+
+         CASE DEFAULT
+            return_value = 1
+            CALL internal_message( 'error', routine_name //                              &
+                                            ': file "' // TRIM( filename ) //            &
+                                            '": file format "' // TRIM( file_format ) // &
+                                            '" not supported' )
+
+      END SELECT
+
+   ENDIF
+
+END FUNCTION dom_write_var
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Find a requested variable 'var_name' and its used dimensions in requested file 'filename'.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE find_var_in_file( filename, var_name, file_format, file_id, var_id, &
+                             is_global, dimensions, return_value )
+
+   CHARACTER(LEN=charlen), INTENT(OUT) ::  file_format  !< file format chosen for file
+   CHARACTER(LEN=*),       INTENT(IN)  ::  filename     !< name of file
+   CHARACTER(LEN=*),       INTENT(IN)  ::  var_name     !< name of variable
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'find_var_in_file'  !< name of routine
+
+   INTEGER(iwp)              ::  d             !< loop index
+   INTEGER(iwp)              ::  dd            !< loop index
+   INTEGER(iwp)              ::  f             !< loop index
+   INTEGER(iwp), INTENT(OUT) ::  file_id       !< file ID
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+   INTEGER(iwp), INTENT(OUT) ::  var_id        !< variable ID
+
+   INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  dim_ids  !< list of dimension IDs used by variable
+
+   LOGICAL              ::  found      !< true if requested variable found in requested file
+   LOGICAL, INTENT(OUT) ::  is_global  !< true if variable is global
+
+   TYPE(dimension_type), DIMENSION(:), ALLOCATABLE, INTENT(OUT) ::  dimensions  !< list of dimensions used by variable
+
+
+   return_value   = 0
+   found = .FALSE.
+
+   DO  f = 1, nf
+      IF ( TRIM( filename ) == TRIM( files(f)%name ) )  THEN
+         file_id     = files(f)%id
+         file_format = files(f)%format
+
+         !-- Search for variable in file
+         DO  d = 1, SIZE( files(f)%variables )
+            IF ( TRIM( var_name ) == TRIM( files(f)%variables(d)%name ) )  THEN
+
+               var_id    = files(f)%variables(d)%id
+               is_global = files(f)%variables(d)%is_global
+
+               ALLOCATE( dim_ids(SIZE( files(f)%variables(d)%dimension_ids )) )
+               ALLOCATE( dimensions(SIZE( files(f)%variables(d)%dimension_ids )) )
+
+               dim_ids = files(f)%variables(d)%dimension_ids
+
+               found = .TRUE.
+               EXIT
+
+            ENDIF
+         ENDDO
+
+         IF ( found )  THEN
+
+            !-- Get list of dimensions used by variable
+            DO  d = 1, SIZE( files(f)%dimensions )
+               DO  dd = 1, SIZE( dim_ids )
+                  IF ( dim_ids(dd) == files(f)%dimensions(d)%id )  THEN
+                     dimensions(dd) = files(f)%dimensions(d)
+                     EXIT
+                  ENDIF
+               ENDDO
+            ENDDO
+
+         ELSE
+
+            !-- If variable was not found, search for a dimension instead
+            DO  d = 1, SIZE( files(f)%dimensions )
+               IF ( TRIM( var_name ) == TRIM( files(f)%dimensions(d)%name ) )  THEN
+
+                  var_id    = files(f)%dimensions(d)%var_id
+                  is_global = .TRUE.
+
+                  ALLOCATE( dimensions(1) )
+
+                  dimensions(1) = files(f)%dimensions(d)
+
+                  found = .TRUE.
+                  EXIT
+
+               ENDIF
+            ENDDO
+
+         ENDIF
+
+         !-- If variable was not found in requested file, return an error
+         IF ( .NOT. found )  THEN
+            return_value = 1
+            CALL internal_message( 'error', routine_name //                       &
+                                            ': variable "' // TRIM( var_name ) // &
+                                            '" not found in file "' // TRIM( filename ) // '"' )
+         ENDIF
+
+         EXIT
+
+      ENDIF  ! file found
+   ENDDO  ! loop over files
+
+   IF ( .NOT. found  .AND.  return_value == 0 )  THEN
+      return_value = 1
+      CALL internal_message( 'error', routine_name //                           &
+                                      ': file "' // TRIM( filename ) //         &
+                                      '" for variable "' // TRIM( var_name ) // &
+                                      '" not found' )
+   ENDIF
+
+END SUBROUTINE find_var_in_file
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Search for masked indices of dimensions within the given bounds ('bounds_start' and
+!> 'bounds_end'). Return the masked indices ('masked_indices') of the dimensions and the bounds of
+!> the masked dimensions containing these indices ('bounds_masked_start' and 'bounds_masked_end').
+!> If no masked indices lie within the given bounds, 'mask_in_bounds' is .FALSE..
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE get_masked_indices_and_masked_dimension_bounds(                                  &
+              dimensions, bounds_start, bounds_end, bounds_masked_start, bounds_masked_end, &
+              masked_indices, mask_in_bounds )
+
+   ! CHARACTER(LEN=*), PARAMETER ::  routine_name = 'get_masked_indices_and_masked_dimension_bounds'  !< name of routine
+
+   INTEGER(iwp) ::  d  !< loop index
+   INTEGER(iwp) ::  i  !< loop index
+   INTEGER(iwp) ::  j  !< count index
+
+   INTEGER(iwp), DIMENSION(:), INTENT(IN)  ::  bounds_end           !< upper bonuds to be searched in
+   INTEGER(iwp), DIMENSION(:), INTENT(OUT) ::  bounds_masked_end    !< upper bounds of masked dimensions within given bounds
+   INTEGER(iwp), DIMENSION(:), INTENT(OUT) ::  bounds_masked_start  !< lower bounds of masked dimensions within given bounds
+   INTEGER(iwp), DIMENSION(:), INTENT(IN)  ::  bounds_start         !< lower bounds to be searched in
+
+   INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE, INTENT(OUT) ::  masked_indices  !< masked indices within given bounds
+
+   LOGICAL, INTENT(OUT) ::  mask_in_bounds  !< true if any masked indices lie within given bounds
+
+   TYPE(dimension_type), DIMENSION(:), INTENT(IN) ::  dimensions  !< dimensions to be searched for masked indices
+
+
+   mask_in_bounds = .TRUE.
+   ALLOCATE( masked_indices(SIZE( dimensions ),0:MAXVAL( dimensions(:)%length - 1 )) )
+   masked_indices = -HUGE( 0_iwp )
+
+   !-- Check for masking and update lower and upper bounds if masked
+   DO  d = 1, SIZE( dimensions )
+
+      IF ( dimensions(d)%is_masked )  THEN
+
+         bounds_masked_end(d) = HUGE( 0_iwp )
+         bounds_masked_start(d) = -HUGE( 0_iwp )
+
+         !-- Find number of masked values within given variable bounds
+         j = 0
+         DO  i = LBOUND( dimensions(d)%masked_indices, DIM=1 ), &
+                 UBOUND( dimensions(d)%masked_indices, DIM=1 )
+
+            !-- Is masked position within given bounds?
+            IF ( dimensions(d)%masked_indices(i) >= bounds_start(d)  .AND.  &
+                 dimensions(d)%masked_indices(i) <= bounds_end(d)           )  THEN
+
+               !-- Save masked position
+               masked_indices(d,i) = dimensions(d)%masked_indices(i)
+               j = j + 1
+
+               !-- Save bounds of mask within given bounds
+               IF ( bounds_masked_start(d) == -HUGE( 0_iwp ) )  bounds_masked_start(d) = i
+               bounds_masked_end(d) = i
+
+            ENDIF
+
+         ENDDO
+
+         !-- Set masked bounds to zero if no masked data lies within bounds
+         IF ( j == 0 )  THEN
+            bounds_masked_start(:) = 0_iwp
+            bounds_masked_end(:) = 0_iwp
+            mask_in_bounds = .FALSE.
+            EXIT
+         ENDIF
+
+      ELSE
+
+         !-- If dimension is not masked, tread all data within bounds as masked
+         bounds_masked_start(d) = bounds_start(d)
+         bounds_masked_end(d) = bounds_end(d)
+         DO  i = bounds_masked_start(d), bounds_masked_end(d)
+            masked_indices(d,i) = i
+         ENDDO
+
+      ENDIF
+
+   ENDDO
+
+END SUBROUTINE get_masked_indices_and_masked_dimension_bounds
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Finalize output.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION dom_finalize_output() RESULT( return_value )
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'dom_finalize_output'  !< name of routine
+
+   INTEGER(iwp) ::  return_value           !< return value
+   INTEGER(iwp) ::  return_value_internal  !< return value from called routines
+   INTEGER(iwp) ::  f                      !< loop index
+
+
+   return_value = 0
+
+   DO  f = 1, nf
+
+      return_value_internal = 0
+
+      SELECT CASE ( TRIM( files(f)%format ) )
+
+         CASE ( 'binary' )
+            CALL binary_finalize( files(f)%id, return_value_internal )
+
+         CASE ( 'netcdf4-serial' )
+            CALL netcdf4_serial_finalize( files(f)%id, return_value_internal )
+
+         CASE ( 'netcdf4-parallel' )
+            CALL netcdf4_parallel_finalize( files(f)%id, return_value_internal )
+
+         CASE DEFAULT
+            return_value_internal = 1
+            CALL internal_message( 'error', routine_name // &
+                                   ': unsupported file format "' // TRIM( files(f)%format ) )
+
+      END SELECT
+
+      IF ( return_value_internal /= 0 )  return_value = return_value_internal
+
+   ENDDO
+
+END FUNCTION dom_finalize_output
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Message routine writing debug information into the debug file
+!> or creating the error message string.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE internal_message( level, string )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  level   !< message importance level
+   CHARACTER(LEN=*), INTENT(IN) ::  string  !< message string
+
+
+   IF ( TRIM( level ) == 'error' )  THEN
+
+      WRITE( internal_error_message, '(A,A)' ) 'DOM ERROR: ', string
+
+   ELSEIF ( TRIM( level ) == 'debug'  .AND.  print_debug_output )  THEN
+
+      WRITE( debug_output_unit, '(A,A)' ) 'DOM DEBUG: ', string
+      FLUSH( debug_output_unit )
+
+   ENDIF
+
+END SUBROUTINE internal_message
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Return the last created error message.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE dom_get_error_message( error_message )
+
+   CHARACTER(LEN=800), INTENT(OUT) ::  error_message         !< return error message to main program
+   CHARACTER(LEN=800)              ::  module_error_message  !< error message created by other module
+
+
+   CALL binary_get_error_message( module_error_message )
+   internal_error_message = TRIM( internal_error_message ) // module_error_message
+
+   CALL netcdf4_serial_get_error_message( module_error_message )
+   internal_error_message = TRIM( internal_error_message ) // module_error_message
+
+   CALL netcdf4_parallel_get_error_message( module_error_message )
+   internal_error_message = TRIM( internal_error_message ) // module_error_message
+
+   error_message = internal_error_message
+
+END SUBROUTINE dom_get_error_message
+
+END MODULE data_output_module
Index: palm/trunk/SOURCE/data_output_netcdf4_parallel_module.f90
===================================================================
--- palm/trunk/SOURCE/data_output_netcdf4_parallel_module.f90	(revision 4070)
+++ palm/trunk/SOURCE/data_output_netcdf4_parallel_module.f90	(revision 4070)
@@ -0,0 +1,787 @@
+!> @file data_output_netcdf4_parallel_module.f90
+!--------------------------------------------------------------------------------------------------!
+! This file is part of the PALM model system.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 2019-2019 Leibniz Universitaet Hannover
+!--------------------------------------------------------------------------------------------------!
+!
+! Current revisions:
+! ------------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
+! Initial revision
+!
+!
+! Authors:
+! --------
+!> @author: Tobias Gronemeier
+!
+! Description:
+! ------------
+!> NetCDF output module to write data to NetCDF files using parallel NetCDF.
+!>
+!> @todo Think of removing 'is_init' as its value can be derived from the return
+!>       value of 'return_value'.
+!--------------------------------------------------------------------------------------------------!
+MODULE data_output_netcdf4_parallel_module
+
+   USE kinds
+
+#if defined( __parallel )
+#if defined( __mpifh )
+   INCLUDE "mpif.h"
+#else
+   USE MPI
+#endif
+#endif
+
+#if defined( __netcdf4_parallel )
+   USE NETCDF
+#endif
+
+   IMPLICIT NONE
+
+   CHARACTER(LEN=800) ::  internal_error_message = ''  !< string containing the last error message
+   CHARACTER(LEN=800) ::  temp_string                  !< dummy string
+
+   INTEGER(iwp) ::  debug_output_unit       !< Fortran Unit Number of the debug-output file
+   INTEGER(iwp) ::  global_id_in_file = -1  !< value of global ID within a file
+
+   LOGICAL ::  print_debug_output = .FALSE.  !< if true, debug output is printed
+
+   SAVE
+
+   PRIVATE
+
+   INTERFACE netcdf4_parallel_init_module
+      MODULE PROCEDURE netcdf4_parallel_init_module
+   END INTERFACE netcdf4_parallel_init_module
+
+   INTERFACE netcdf4_parallel_open_file
+      MODULE PROCEDURE netcdf4_parallel_open_file
+   END INTERFACE netcdf4_parallel_open_file
+
+   INTERFACE netcdf4_parallel_init_dimension
+      MODULE PROCEDURE netcdf4_parallel_init_dimension
+   END INTERFACE netcdf4_parallel_init_dimension
+
+   INTERFACE netcdf4_parallel_init_variable
+      MODULE PROCEDURE netcdf4_parallel_init_variable
+   END INTERFACE netcdf4_parallel_init_variable
+
+   INTERFACE netcdf4_parallel_write_attribute
+      MODULE PROCEDURE netcdf4_parallel_write_attribute
+   END INTERFACE netcdf4_parallel_write_attribute
+
+   INTERFACE netcdf4_parallel_init_end
+      MODULE PROCEDURE netcdf4_parallel_init_end
+   END INTERFACE netcdf4_parallel_init_end
+
+   INTERFACE netcdf4_parallel_write_variable
+      MODULE PROCEDURE netcdf4_parallel_write_variable
+   END INTERFACE netcdf4_parallel_write_variable
+
+   INTERFACE netcdf4_parallel_finalize
+      MODULE PROCEDURE netcdf4_parallel_finalize
+   END INTERFACE netcdf4_parallel_finalize
+
+   INTERFACE netcdf4_parallel_get_error_message
+      MODULE PROCEDURE netcdf4_parallel_get_error_message
+   END INTERFACE netcdf4_parallel_get_error_message
+
+   PUBLIC &
+      netcdf4_parallel_finalize, &
+      netcdf4_parallel_get_error_message, &
+      netcdf4_parallel_init_dimension, &
+      netcdf4_parallel_init_end, &
+      netcdf4_parallel_init_module, &
+      netcdf4_parallel_init_variable, &
+      netcdf4_parallel_open_file, &
+      netcdf4_parallel_write_attribute, &
+      netcdf4_parallel_write_variable
+
+
+CONTAINS
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize data-output module.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_parallel_init_module( program_debug_output_unit, debug_output, dom_global_id )
+
+   INTEGER(iwp), INTENT(IN) ::  dom_global_id              !< global id within a file defined by DOM
+   INTEGER(iwp), INTENT(IN) ::  program_debug_output_unit  !< file unit number for debug output
+
+   LOGICAL, INTENT(IN) ::  debug_output  !< if true, debug output is printed
+
+
+   debug_output_unit = program_debug_output_unit
+
+   print_debug_output = debug_output
+
+   global_id_in_file = dom_global_id
+
+
+END SUBROUTINE netcdf4_parallel_init_module
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Open netcdf file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_parallel_open_file( filename, file_id, return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  filename  !< name of file
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_parallel_open_file'  !< name of this routine
+
+   INTEGER(iwp), INTENT(OUT) ::  file_id       !< file ID
+   INTEGER(iwp)              ::  nc_stat       !< netcdf return value
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+
+#if defined( __parallel ) && defined( __netcdf4_parallel )
+   return_value = 0
+
+   !-- Open new file
+   CALL internal_message( 'debug', routine_name // ': create file "' // TRIM( filename ) // '"' )
+
+   nc_stat = NF90_CREATE( TRIM( filename ),                                       &
+                          IOR( NF90_NOCLOBBER, IOR( NF90_NETCDF4, NF90_MPIIO ) ), &
+                          file_id, COMM = MPI_COMM_WORLD, INFO = MPI_INFO_NULL )
+
+   IF ( nc_stat /= NF90_NOERR )  THEN
+      return_value = 1
+      CALL internal_message( 'error', routine_name // ': NetCDF error while opening file "' // &
+                                      TRIM( filename ) // '": ' // NF90_STRERROR( nc_stat ) )
+   ENDIF
+#else
+   file_id = -1
+   nc_stat = 0
+   return_value = 1
+   CALL internal_message( 'error', routine_name //                      &
+                          ': pre-processor directives "__parallel" ' // &
+                          'and/or "__netcdf4_parallel" not given. ' //  &
+                          'Using parallel NetCDF4 output not possible' )
+#endif
+
+END SUBROUTINE netcdf4_parallel_open_file
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Write attribute to netcdf file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_parallel_write_attribute( file_id, var_id, att_name, att_value_char, &
+                 att_value_int8, att_value_int16, att_value_int32,                      &
+                 att_value_real32, att_value_real64, return_value )
+
+   CHARACTER(LEN=*), INTENT(IN)           ::  att_name        !< name of attribute
+   CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  att_value_char  !< value of attribute
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_parallel_write_attribute'  !< name of this routine
+
+   INTEGER(iwp) ::  nc_stat    !< netcdf return value
+   INTEGER(iwp) ::  target_id  !< ID of target which gets attribute (either global or var_id)
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+   INTEGER(iwp), INTENT(IN)  ::  var_id        !< variable ID
+
+   INTEGER(KIND=1), INTENT(IN), OPTIONAL ::  att_value_int8   !< value of attribute
+   INTEGER(KIND=2), INTENT(IN), OPTIONAL ::  att_value_int16  !< value of attribute
+   INTEGER(KIND=4), INTENT(IN), OPTIONAL ::  att_value_int32  !< value of attribute
+
+   REAL(KIND=4), INTENT(IN), OPTIONAL ::  att_value_real32  !< value of attribute
+   REAL(KIND=8), INTENT(IN), OPTIONAL ::  att_value_real64  !< value of attribute
+
+
+#if defined( __netcdf4_parallel )
+   return_value = 0
+
+   IF ( var_id == global_id_in_file )  THEN
+      target_id = NF90_GLOBAL
+   ELSE
+      target_id = var_id
+   ENDIF
+
+   CALL internal_message( 'debug', &
+                          routine_name // ': write attribute "' // TRIM( att_name ) // '"' )
+
+   IF ( PRESENT( att_value_char ) )  THEN
+      nc_stat = NF90_PUT_ATT( file_id, target_id, TRIM( att_name ), TRIM( att_value_char ) )
+   ELSEIF ( PRESENT( att_value_int8 ) )  THEN
+      nc_stat = NF90_PUT_ATT( file_id, target_id, TRIM( att_name ), att_value_int8 )
+   ELSEIF ( PRESENT( att_value_int16 ) )  THEN
+      nc_stat = NF90_PUT_ATT( file_id, target_id, TRIM( att_name ), att_value_int16 )
+   ELSEIF ( PRESENT( att_value_int32 ) )  THEN
+      nc_stat = NF90_PUT_ATT( file_id, target_id, TRIM( att_name ), att_value_int32 )
+   ELSEIF ( PRESENT( att_value_real32 ) )  THEN
+      nc_stat = NF90_PUT_ATT( file_id, target_id, TRIM( att_name ), att_value_real32 )
+   ELSEIF ( PRESENT( att_value_real64 ) )  THEN
+      nc_stat = NF90_PUT_ATT( file_id, target_id, TRIM( att_name ), att_value_real64 )
+   ELSE
+      return_value = 1
+      CALL internal_message( 'error', TRIM( routine_name ) // &
+                                      ': attribute "' // TRIM( att_name ) // '": no value given' )
+   ENDIF
+
+   IF ( return_value == 0 )  THEN
+      IF ( nc_stat /= NF90_NOERR )  THEN
+         return_value = 1
+         CALL internal_message( 'error',                                       &
+                 routine_name // ': NetCDF error while writing attribute "' // &
+                 TRIM( att_name ) // '": ' // NF90_STRERROR( nc_stat ) )
+      ENDIF
+   ENDIF
+#else
+   return_value = 1
+#endif
+
+END SUBROUTINE netcdf4_parallel_write_attribute
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize dimension.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_parallel_init_dimension( file_id, dim_id, var_id, &
+              dim_name, dim_type, dim_length, is_init, return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  dim_name  !< name of dimension
+   CHARACTER(LEN=*), INTENT(IN) ::  dim_type  !< data type of dimension
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_parallel_init_dimension'  !< name of this routine
+
+   INTEGER(iwp), INTENT(OUT) ::  dim_id         !< dimension ID
+   INTEGER(iwp), INTENT(IN)  ::  dim_length     !< length of dimension
+   INTEGER(iwp), INTENT(IN)  ::  file_id        !< file ID
+   INTEGER(iwp)              ::  nc_dim_length  !< length of dimension
+   INTEGER(iwp)              ::  nc_stat        !< netcdf return value
+   INTEGER(iwp), INTENT(OUT) ::  return_value   !< return value
+   INTEGER(iwp), INTENT(OUT) ::  var_id         !< variable ID
+
+   LOGICAL, INTENT(OUT) ::  is_init  !< true if dimension is initialized
+
+
+#if defined( __netcdf4_parallel )
+   return_value = 0
+   var_id = -1
+
+   CALL internal_message( 'debug', routine_name // ': init dimension "' // TRIM( dim_name ) // '"' )
+
+   !-- Check if dimension is unlimited
+   IF ( dim_length < 0 )  THEN
+      nc_dim_length = NF90_UNLIMITED
+   ELSE
+      nc_dim_length = dim_length
+   ENDIF
+
+   !-- Define dimension in file
+   nc_stat = NF90_DEF_DIM( file_id, dim_name, nc_dim_length, dim_id )
+
+   IF ( nc_stat == NF90_NOERR )  THEN
+
+      !-- Define variable holding dimension values in file
+      CALL netcdf4_parallel_init_variable( file_id, var_id, dim_name, dim_type, (/dim_id/), &
+                                           is_init, is_global=.TRUE., return_value=return_value )
+
+   ELSE
+      return_value = 1
+      is_init = .FALSE.
+      CALL internal_message( 'error', routine_name //                                    &
+                                      ': NetCDF error while initializing dimension "' // &
+                                      TRIM( dim_name ) // '": ' // NF90_STRERROR( nc_stat ) )
+   ENDIF
+#else
+   return_value = 1
+   var_id = -1
+   dim_id = -1
+   is_init = .FALSE.
+#endif
+
+END SUBROUTINE netcdf4_parallel_init_dimension
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize variable.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_parallel_init_variable( file_id, var_id, var_name, var_type, var_dim_ids, &
+                                           is_init, is_global, return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  var_name  !< name of variable
+   CHARACTER(LEN=*), INTENT(IN) ::  var_type  !< data type of variable
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_parallel_init_variable'  !< name of this routine
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp)              ::  nc_stat       !< netcdf return value
+   INTEGER(iwp)              ::  nc_var_type   !< netcdf data type
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+   INTEGER(iwp), INTENT(OUT) ::  var_id        !< variable ID
+
+   INTEGER(iwp), DIMENSION(:), INTENT(IN) ::  var_dim_ids  !< list of dimension IDs used by variable
+
+   LOGICAL, INTENT(IN)  ::  is_global  !< true if variable is global (same on all PE)
+   LOGICAL, INTENT(OUT) ::  is_init    !< true if variable is initialized
+
+
+#if defined( __netcdf4_parallel )
+   return_value = 0
+
+   WRITE( temp_string, * ) is_global
+   CALL internal_message( 'debug', routine_name // ': init variable "' // TRIM( var_name ) // &
+                                   '" ( is_global = ' // TRIM( temp_string ) // ')' )
+
+   nc_var_type = get_netcdf_data_type( var_type )
+
+   IF ( nc_var_type /= -1_iwp )  THEN
+
+      !-- Define variable in file
+      nc_stat = NF90_DEF_VAR( file_id, var_name, nc_var_type, var_dim_ids, var_id )
+
+      IF ( nc_stat == NF90_NOERR )  THEN
+         IF ( is_global )  THEN
+            nc_stat = NF90_VAR_PAR_ACCESS( file_id, var_id, NF90_INDEPENDENT )
+         ELSE
+            nc_stat = NF90_VAR_PAR_ACCESS( file_id, var_id, NF90_COLLECTIVE )
+         ENDIF
+      ENDIF
+
+      IF ( nc_stat /= NF90_NOERR)  THEN
+         return_value = 1
+         CALL internal_message( 'error', routine_name //                                   &
+                                         ': NetCDF error while initializing variable "' // &
+                                         TRIM( var_name ) // '": ' // NF90_STRERROR( nc_stat ) )
+      ENDIF
+
+   ELSE
+      return_value = 1
+   ENDIF
+
+   is_init = return_value == 0
+#else
+   return_value = 1
+   var_id = -1
+   is_init = .FALSE.
+#endif
+
+END SUBROUTINE netcdf4_parallel_init_variable
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Leave file definition state.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_parallel_init_end( file_id, return_value )
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_parallel_init_end'  !< name of this routine
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp)              ::  nc_stat       !< netcdf return value
+   INTEGER(iwp)              ::  old_mode      !< previous netcdf fill mode
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+
+#if defined( __netcdf4_parallel )
+   return_value = 0
+
+   WRITE( temp_string, * ) file_id
+   CALL internal_message( 'debug',        &
+                          routine_name // &
+                          ': finalize file definition (file_id=' // TRIM( temp_string ) // ')' )
+
+   !-- Set general no fill, otherwise the performance drops significantly
+   nc_stat = NF90_SET_FILL( file_id, NF90_NOFILL, old_mode )
+
+   IF ( nc_stat == NF90_NOERR )  THEN
+      nc_stat = NF90_ENDDEF( file_id )
+   ENDIF
+
+   IF ( nc_stat /= NF90_NOERR )  THEN
+      return_value = 1
+      CALL internal_message( 'error', routine_name // ': NetCDF error: ' // &
+                                      NF90_STRERROR( nc_stat ) )
+   ENDIF
+#else
+   return_value = 1
+#endif
+
+END SUBROUTINE netcdf4_parallel_init_end
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Write variable of different kind into netcdf file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_parallel_write_variable(                               &
+              file_id, var_id, bounds_start, bounds_end, bounds_origin,   &
+              do_output, is_global,                                       &
+              var_int8_0d,   var_int8_1d,   var_int8_2d,   var_int8_3d,   &
+              var_int16_0d,  var_int16_1d,  var_int16_2d,  var_int16_3d,  &
+              var_int32_0d,  var_int32_1d,  var_int32_2d,  var_int32_3d,  &
+              var_intwp_0d,  var_intwp_1d,  var_intwp_2d,  var_intwp_3d,  &
+              var_real32_0d, var_real32_1d, var_real32_2d, var_real32_3d, &
+              var_real64_0d, var_real64_1d, var_real64_2d, var_real64_3d, &
+              var_realwp_0d, var_realwp_1d, var_realwp_2d, var_realwp_3d, &
+              return_value )
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_parallel_write_variable'  !< name of this routine
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp)              ::  myid = 0      !< id number of processor element
+   INTEGER(iwp)              ::  nc_stat       !< netcdf return value
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+   INTEGER(iwp), INTENT(IN)  ::  var_id        !< variable ID
+
+   INTEGER(iwp), DIMENSION(:), INTENT(IN)  ::  bounds_origin  !< starting index of each dimension
+   INTEGER(iwp), DIMENSION(:), INTENT(IN)  ::  bounds_end     !< ending index of variable
+   INTEGER(iwp), DIMENSION(:), INTENT(IN)  ::  bounds_start   !< starting index of variable
+   INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  value_count    !< count of values along each dimension to be written
+
+   INTEGER(KIND=1), POINTER, INTENT(IN), OPTIONAL                               ::  var_int8_0d  !< output variable
+   INTEGER(KIND=1), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_int8_1d  !< output variable
+   INTEGER(KIND=1), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_int8_2d  !< output variable
+   INTEGER(KIND=1), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_int8_3d  !< output variable
+
+   INTEGER(KIND=2), POINTER, INTENT(IN), OPTIONAL                               ::  var_int16_0d  !< output variable
+   INTEGER(KIND=2), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_int16_1d  !< output variable
+   INTEGER(KIND=2), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_int16_2d  !< output variable
+   INTEGER(KIND=2), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_int16_3d  !< output variable
+
+   INTEGER(KIND=4), POINTER, INTENT(IN), OPTIONAL                               ::  var_int32_0d  !< output variable
+   INTEGER(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_int32_1d  !< output variable
+   INTEGER(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_int32_2d  !< output variable
+   INTEGER(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_int32_3d  !< output variable
+
+   INTEGER(iwp), POINTER, INTENT(IN), OPTIONAL                               ::  var_intwp_0d  !< output variable
+   INTEGER(iwp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_intwp_1d  !< output variable
+   INTEGER(iwp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_intwp_2d  !< output variable
+   INTEGER(iwp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_intwp_3d  !< output variable
+
+   LOGICAL, INTENT(IN) ::  do_output  !< if false, set count to 0 and do no output
+   LOGICAL, INTENT(IN) ::  is_global  !< true if variable is global (same on all PE)
+
+   REAL(KIND=4), POINTER, INTENT(IN), OPTIONAL                               ::  var_real32_0d  !< output variable
+   REAL(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_real32_1d  !< output variable
+   REAL(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_real32_2d  !< output variable
+   REAL(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_real32_3d  !< output variable
+
+   REAL(KIND=8), POINTER, INTENT(IN), OPTIONAL                               ::  var_real64_0d  !< output variable
+   REAL(KIND=8), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_real64_1d  !< output variable
+   REAL(KIND=8), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_real64_2d  !< output variable
+   REAL(KIND=8), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_real64_3d  !< output variable
+
+   REAL(wp), POINTER, INTENT(IN), OPTIONAL                               ::  var_realwp_0d  !< output variable
+   REAL(wp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_realwp_1d  !< output variable
+   REAL(wp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_realwp_2d  !< output variable
+   REAL(wp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_realwp_3d  !< output variable
+
+
+#if defined( __netcdf4_parallel )
+   return_value = 0
+
+#if defined( __parallel )
+   CALL MPI_COMM_RANK( MPI_COMM_WORLD, myid, return_value )
+   IF ( return_value /= 0 )  THEN
+      CALL internal_message( 'error', routine_name // ': MPI_COMM_RANK error' )
+   ENDIF
+#endif
+
+   IF ( return_value == 0  .AND.  ( .NOT. is_global  .OR.  myid == 0 ) )  THEN
+
+      WRITE( temp_string, * ) var_id
+      CALL internal_message( 'debug', routine_name // ': write variable ' // TRIM( temp_string ) )
+
+      ALLOCATE( value_count(SIZE( bounds_start )) )
+
+      IF ( do_output ) THEN
+         value_count = bounds_end - bounds_start + 1
+      ELSE
+         value_count = 0
+      END IF
+
+      !-- 8bit integer output
+      IF ( PRESENT( var_int8_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_int8_0d /),       &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int8_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int8_1d,             &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int8_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int8_2d,             &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int8_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int8_3d,             &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      !-- 16bit integer output
+      ELSEIF ( PRESENT( var_int16_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_int16_0d /),      &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int16_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int16_1d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int16_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int16_2d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int16_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int16_3d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      !-- 32bit integer output
+      ELSEIF ( PRESENT( var_int32_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_int32_0d /),      &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int32_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int32_1d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int32_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int32_2d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int32_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int32_3d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      !-- working-precision integer output
+      ELSEIF ( PRESENT( var_intwp_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_intwp_0d /),      &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_intwp_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_intwp_1d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_intwp_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_intwp_2d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_intwp_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_intwp_3d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      !-- 32bit real output
+      ELSEIF ( PRESENT( var_real32_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_real32_0d /),     &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_real32_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_real32_1d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_real32_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_real32_2d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_real32_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_real32_3d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      !-- 64bit real output
+      ELSEIF ( PRESENT( var_real64_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_real64_0d /),     &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_real64_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_real64_1d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_real64_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_real64_2d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_real64_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_real64_3d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      !-- working-precision real output
+      ELSEIF ( PRESENT( var_realwp_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_realwp_0d /),     &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_realwp_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_realwp_1d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_realwp_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_realwp_2d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_realwp_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_realwp_3d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSE
+         return_value = 1
+         nc_stat = NF90_NOERR
+         WRITE( temp_string, '(": variable_id=",I6,", file_id=",I6)' ) var_id, file_id
+         CALL internal_message( 'error', routine_name // &
+                                TRIM( temp_string ) //   &
+                                ': no output values given' )
+      ENDIF
+
+      !-- Check for errors
+      IF ( nc_stat /= NF90_NOERR )  THEN
+         return_value = 1
+         WRITE( temp_string, * ) 'variable_id=', var_id, '; file_id=', file_id, &
+                                 ', bounds_start=', bounds_start,               &
+                                 ', bounds_end=', bounds_end,                   &
+                                 ', count=', value_count
+         CALL internal_message( 'error', routine_name //                                    &
+                                         ': error while writing ' // TRIM( temp_string ) // &
+                                         ': NetCDF error: ' // NF90_STRERROR( nc_stat ) )
+      ENDIF
+
+   ENDIF
+#else
+   return_value = 1
+#endif
+
+END SUBROUTINE netcdf4_parallel_write_variable
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Close netcdf file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_parallel_finalize( file_id, return_value )
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_parallel_finalize'  !< name of routine
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp)              ::  nc_stat       !< netcdf return value
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+
+#if defined( __netcdf4_parallel )
+   WRITE( temp_string, * ) file_id
+   CALL internal_message( 'debug', routine_name // &
+                                   ': close file (file_id=' // TRIM( temp_string ) // ')' )
+
+   nc_stat = NF90_CLOSE( file_id )
+   IF ( nc_stat == NF90_NOERR )  THEN
+      return_value = 0
+   ELSE
+      return_value = 1
+      CALL internal_message( 'error', routine_name // &
+                                      ': NetCDF error: ' // NF90_STRERROR( nc_stat ) )
+   ENDIF
+#else
+   return_value = 1
+#endif
+
+END SUBROUTINE netcdf4_parallel_finalize
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Convert data_type string into netcdf data type value.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION get_netcdf_data_type( data_type ) RESULT( return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  data_type  !< requested data type
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'get_netcdf_data_type'  !< name of this routine
+
+   INTEGER(iwp) ::  return_value  !< netcdf data type
+
+
+   SELECT CASE ( TRIM( data_type ) )
+
+#if defined( __netcdf4_parallel )
+      CASE ( 'char' )
+         return_value = NF90_CHAR
+
+      CASE ( 'int8' )
+         return_value = NF90_BYTE
+
+      CASE ( 'int16' )
+         return_value = NF90_SHORT
+
+      CASE ( 'int32' )
+         return_value = NF90_INT
+
+      CASE ( 'real32' )
+         return_value = NF90_FLOAT
+
+      CASE ( 'real64' )
+         return_value = NF90_DOUBLE
+#endif
+
+      CASE DEFAULT
+         CALL internal_message( 'error', routine_name // &
+                                         ': data type unknown (' // TRIM( data_type ) // ')' )
+         return_value = -1_iwp
+
+   END SELECT
+
+END FUNCTION get_netcdf_data_type
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Message routine writing debug information into the debug file
+!> or creating the error message string.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE internal_message( level, string )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  level   !< message importance level
+   CHARACTER(LEN=*), INTENT(IN) ::  string  !< message string
+
+
+   IF ( TRIM( level ) == 'error' )  THEN
+
+      WRITE( internal_error_message, '(A,A)' ) 'DOM ERROR: ', string
+
+   ELSEIF ( TRIM( level ) == 'debug'  .AND.  print_debug_output )  THEN
+
+      WRITE( debug_output_unit, '(A,A)' ) 'DOM DEBUG: ', string
+      FLUSH( debug_output_unit )
+
+   ENDIF
+
+END SUBROUTINE internal_message
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Return the last created error message.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_parallel_get_error_message( error_message )
+
+   CHARACTER(LEN=800), INTENT(OUT) ::  error_message  !< return error message to main program
+
+
+   error_message = internal_error_message
+
+END SUBROUTINE netcdf4_parallel_get_error_message
+
+
+END MODULE data_output_netcdf4_parallel_module
Index: palm/trunk/SOURCE/data_output_netcdf4_serial_module.f90
===================================================================
--- palm/trunk/SOURCE/data_output_netcdf4_serial_module.f90	(revision 4070)
+++ palm/trunk/SOURCE/data_output_netcdf4_serial_module.f90	(revision 4070)
@@ -0,0 +1,780 @@
+!> @file data_output_netcdf4_serial_module.f90
+!--------------------------------------------------------------------------------------------------!
+! This file is part of the PALM model system.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 2019-2019 Leibniz Universitaet Hannover
+!--------------------------------------------------------------------------------------------------!
+!
+! Current revisions:
+! ------------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
+! Initial revision
+!
+!
+! Authors:
+! --------
+!> @author: Tobias Gronemeier
+!
+! Description:
+! ------------
+!> Netcdf output module to write data to NetCDF files by single processor
+!>
+!> @todo Think of removing 'is_init' as its value can be derived from the return
+!>       value of 'return_value'.
+!--------------------------------------------------------------------------------------------------!
+MODULE data_output_netcdf4_serial_module
+
+   USE kinds
+
+#if defined( __parallel )
+#if defined( __mpifh )
+   INCLUDE "mpif.h"
+#else
+   USE MPI
+#endif
+#endif
+
+#if defined( __netcdf4 )
+   USE NETCDF
+#endif
+
+   IMPLICIT NONE
+
+   CHARACTER(LEN=800) ::  internal_error_message = ''  !< string containing the last error message
+   CHARACTER(LEN=800) ::  temp_string                  !< dummy string
+
+   INTEGER(iwp) ::  debug_output_unit       !< Fortran Unit Number of the debug-output file
+   INTEGER(iwp) ::  global_id_in_file = -1  !< value of global ID within a file
+
+   LOGICAL ::  print_debug_output = .FALSE.  !< if true, debug output is printed
+
+   SAVE
+
+   PRIVATE
+
+   INTERFACE netcdf4_serial_init_module
+      MODULE PROCEDURE netcdf4_serial_init_module
+   END INTERFACE netcdf4_serial_init_module
+
+   INTERFACE netcdf4_serial_open_file
+      MODULE PROCEDURE netcdf4_serial_open_file
+   END INTERFACE netcdf4_serial_open_file
+
+   INTERFACE netcdf4_serial_init_dimension
+      MODULE PROCEDURE netcdf4_serial_init_dimension
+   END INTERFACE netcdf4_serial_init_dimension
+
+   INTERFACE netcdf4_serial_init_variable
+      MODULE PROCEDURE netcdf4_serial_init_variable
+   END INTERFACE netcdf4_serial_init_variable
+
+   INTERFACE netcdf4_serial_write_attribute
+      MODULE PROCEDURE netcdf4_serial_write_attribute
+   END INTERFACE netcdf4_serial_write_attribute
+
+   INTERFACE netcdf4_serial_init_end
+      MODULE PROCEDURE netcdf4_serial_init_end
+   END INTERFACE netcdf4_serial_init_end
+
+   INTERFACE netcdf4_serial_write_variable
+      MODULE PROCEDURE netcdf4_serial_write_variable
+   END INTERFACE netcdf4_serial_write_variable
+
+   INTERFACE netcdf4_serial_finalize
+      MODULE PROCEDURE netcdf4_serial_finalize
+   END INTERFACE netcdf4_serial_finalize
+
+   INTERFACE netcdf4_serial_get_error_message
+      MODULE PROCEDURE netcdf4_serial_get_error_message
+   END INTERFACE netcdf4_serial_get_error_message
+
+   PUBLIC &
+   netcdf4_serial_finalize, &
+   netcdf4_serial_get_error_message, &
+   netcdf4_serial_init_dimension, &
+   netcdf4_serial_init_end, &
+   netcdf4_serial_init_module, &
+   netcdf4_serial_init_variable, &
+   netcdf4_serial_open_file, &
+   netcdf4_serial_write_attribute, &
+   netcdf4_serial_write_variable
+
+
+CONTAINS
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize data-output module.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_serial_init_module( program_debug_output_unit, debug_output, dom_global_id )
+
+   INTEGER(iwp), INTENT(IN) ::  dom_global_id              !< global id within a file defined by DOM
+   INTEGER(iwp), INTENT(IN) ::  program_debug_output_unit  !< file unit number for debug output
+
+   LOGICAL, INTENT(IN) ::  debug_output  !< if true, debug output is printed
+
+
+   debug_output_unit = program_debug_output_unit
+
+   print_debug_output = debug_output
+
+   global_id_in_file = dom_global_id
+
+END SUBROUTINE netcdf4_serial_init_module
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Open netcdf file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_serial_open_file( filename, file_id, return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  filename  !< name of file
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_serial_open_file'  !< name of this routine
+
+   INTEGER(iwp), INTENT(OUT) ::  file_id       !< file ID
+   INTEGER(iwp)              ::  nc_stat       !< netcdf return value
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+
+#if defined( __netcdf4 )
+   return_value = 0
+
+   !-- Open new file
+   CALL internal_message( 'debug', routine_name // ': create file "' // TRIM( filename ) // '"' )
+
+   nc_stat = NF90_CREATE( TRIM( filename ),                    &
+                          IOR( NF90_NOCLOBBER, NF90_NETCDF4 ), &
+                          file_id )
+
+   IF ( nc_stat /= NF90_NOERR )  THEN
+      return_value = 1
+      CALL internal_message( 'error', routine_name // ': NetCDF error while opening file "' // &
+                                      TRIM( filename ) // '": ' // NF90_STRERROR( nc_stat ) )
+   ENDIF
+#else
+   file_id = -1
+   return_value = 1
+   CALL internal_message( 'error', routine_name //                               &
+                          ': pre-processor directive "__netcdf4" not given. ' // &
+                          'Using NetCDF4 output not possible' )
+#endif
+
+END SUBROUTINE netcdf4_serial_open_file
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Write attribute to netcdf file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_serial_write_attribute( file_id, var_id, att_name, att_value_char, &
+                 att_value_int8, att_value_int16, att_value_int32,                    &
+                 att_value_real32, att_value_real64, return_value )
+
+   CHARACTER(LEN=*), INTENT(IN)           ::  att_name        !< name of attribute
+   CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  att_value_char  !< value of attribute
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_serial_write_attribute'  !< name of this routine
+
+   INTEGER(iwp) ::  nc_stat    !< netcdf return value
+   INTEGER(iwp) ::  target_id  !< ID of target which gets attribute (either global or var_id)
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+   INTEGER(iwp), INTENT(IN)  ::  var_id        !< variable ID
+
+   INTEGER(KIND=1), INTENT(IN), OPTIONAL ::  att_value_int8   !< value of attribute
+   INTEGER(KIND=2), INTENT(IN), OPTIONAL ::  att_value_int16  !< value of attribute
+   INTEGER(KIND=4), INTENT(IN), OPTIONAL ::  att_value_int32  !< value of attribute
+
+   REAL(KIND=4), INTENT(IN), OPTIONAL ::  att_value_real32  !< value of attribute
+   REAL(KIND=8), INTENT(IN), OPTIONAL ::  att_value_real64  !< value of attribute
+
+
+#if defined( __netcdf4 )
+   return_value = 0
+
+   IF ( var_id == global_id_in_file )  THEN
+      target_id = NF90_GLOBAL
+   ELSE
+      target_id = var_id
+   ENDIF
+
+   CALL internal_message( 'debug', &
+                          routine_name // ': write attribute "' // TRIM( att_name ) // '"' )
+
+   IF ( PRESENT( att_value_char ) )  THEN
+      nc_stat = NF90_PUT_ATT( file_id, target_id, TRIM( att_name ), TRIM( att_value_char ) )
+   ELSEIF ( PRESENT( att_value_int8 ) )  THEN
+      nc_stat = NF90_PUT_ATT( file_id, target_id, TRIM( att_name ), att_value_int8 )
+   ELSEIF ( PRESENT( att_value_int16 ) )  THEN
+      nc_stat = NF90_PUT_ATT( file_id, target_id, TRIM( att_name ), att_value_int16 )
+   ELSEIF ( PRESENT( att_value_int32 ) )  THEN
+      nc_stat = NF90_PUT_ATT( file_id, target_id, TRIM( att_name ), att_value_int32 )
+   ELSEIF ( PRESENT( att_value_real32 ) )  THEN
+      nc_stat = NF90_PUT_ATT( file_id, target_id, TRIM( att_name ), att_value_real32 )
+   ELSEIF ( PRESENT( att_value_real64 ) )  THEN
+      nc_stat = NF90_PUT_ATT( file_id, target_id, TRIM( att_name ), att_value_real64 )
+   ELSE
+      return_value = 1
+      CALL internal_message( 'error', TRIM( routine_name ) // &
+                                      ': attribute "' // TRIM( att_name ) // '": no value given' )
+   ENDIF
+
+   IF ( return_value == 0 )  THEN
+      IF ( nc_stat /= NF90_NOERR )  THEN
+         return_value = 1
+         CALL internal_message( 'error',                                       &
+                 routine_name // ': NetCDF error while writing attribute "' // &
+                 TRIM( att_name ) // '": ' // NF90_STRERROR( nc_stat ) )
+      ENDIF
+   ENDIF
+#else
+   return_value = 1
+#endif
+
+END SUBROUTINE netcdf4_serial_write_attribute
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize dimension.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_serial_init_dimension( file_id, dim_id, var_id, &
+              dim_name, dim_type, dim_length, is_init, return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  dim_name  !< name of dimension
+   CHARACTER(LEN=*), INTENT(IN) ::  dim_type  !< data type of dimension
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_serial_init_dimension'  !< name of this routine
+
+   INTEGER(iwp), INTENT(OUT) ::  dim_id         !< dimension ID
+   INTEGER(iwp), INTENT(IN)  ::  dim_length     !< length of dimension
+   INTEGER(iwp), INTENT(IN)  ::  file_id        !< file ID
+   INTEGER(iwp)              ::  nc_dim_length  !< length of dimension
+   INTEGER(iwp)              ::  nc_stat        !< netcdf return value
+   INTEGER(iwp), INTENT(OUT) ::  return_value   !< return value
+   INTEGER(iwp), INTENT(OUT) ::  var_id         !< variable ID
+
+   LOGICAL, INTENT(OUT) ::  is_init  !< true if dimension is initialized
+
+
+#if defined( __netcdf4 )
+   return_value = 0
+   var_id = -1
+
+   CALL internal_message( 'debug', routine_name // ': init dimension "' // TRIM( dim_name ) // '"' )
+
+   !-- Check if dimension is unlimited
+   IF ( dim_length < 0 )  THEN
+      nc_dim_length = NF90_UNLIMITED
+   ELSE
+      nc_dim_length = dim_length
+   ENDIF
+
+   !-- Define dimension in file
+   nc_stat = NF90_DEF_DIM( file_id, dim_name, nc_dim_length, dim_id )
+
+   IF ( nc_stat == NF90_NOERR )  THEN
+
+      !-- Define variable holding dimension values in file
+      CALL netcdf4_serial_init_variable( file_id, var_id, dim_name, dim_type, (/dim_id/), &
+                                         is_init, is_global=.TRUE., return_value=return_value )
+
+   ELSE
+      return_value = 1
+      is_init = .FALSE.
+      CALL internal_message( 'error', routine_name //                                    &
+                                      ': NetCDF error while initializing dimension "' // &
+                                      TRIM( dim_name ) // '": ' // NF90_STRERROR( nc_stat ) )
+   ENDIF
+#else
+   return_value = 1
+   var_id = -1
+   dim_id = -1
+   is_init = .FALSE.
+#endif
+
+END SUBROUTINE netcdf4_serial_init_dimension
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialize variable.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_serial_init_variable( file_id, var_id, var_name, var_type, var_dim_ids, &
+                                         is_init, is_global, return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  var_name  !< name of variable
+   CHARACTER(LEN=*), INTENT(IN) ::  var_type  !< data type of variable
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_serial_init_variable'  !< name of this routine
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp)              ::  nc_stat       !< netcdf return value
+   INTEGER(iwp)              ::  nc_var_type   !< netcdf data type
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+   INTEGER(iwp), INTENT(OUT) ::  var_id        !< variable ID
+
+   INTEGER(iwp), DIMENSION(:), INTENT(IN) ::  var_dim_ids  !< list of dimension IDs used by variable
+
+   LOGICAL, INTENT(IN)  ::  is_global  !< true if variable is global (same on all PE)
+   LOGICAL, INTENT(OUT) ::  is_init    !< true if variable is initialized
+
+
+#if defined( __netcdf4 )
+   return_value = 0
+
+   WRITE( temp_string, * ) is_global
+   CALL internal_message( 'debug', routine_name // ': init variable "' // TRIM( var_name ) // &
+                                   '" ( is_global = ' // TRIM( temp_string ) // ')' )
+
+   nc_var_type = get_netcdf_data_type( var_type )
+
+   IF ( nc_var_type /= -1_iwp )  THEN
+
+      !-- Define variable in file
+      nc_stat = NF90_DEF_VAR( file_id, var_name, nc_var_type, var_dim_ids, var_id )
+
+      IF ( is_global )  CONTINUE  ! Added to be consistent to parallel output
+
+      IF ( nc_stat /= NF90_NOERR )  THEN
+         return_value = 1
+         CALL internal_message( 'error', routine_name //                                   &
+                                         ': NetCDF error while initializing variable "' // &
+                                         TRIM( var_name ) // '": ' // NF90_STRERROR( nc_stat ) )
+      ENDIF
+
+   ELSE
+      return_value = 1
+   ENDIF
+
+   is_init = return_value == 0
+#else
+   return_value = 1
+   var_id = -1
+   is_init = .FALSE.
+#endif
+
+END SUBROUTINE netcdf4_serial_init_variable
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Leave file definition state.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_serial_init_end( file_id, return_value )
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_serial_init_end'  !< name of this routine
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp)              ::  nc_stat       !< netcdf return value
+   INTEGER(iwp)              ::  old_mode      !< previous netcdf fill mode
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+
+#if defined( __netcdf4 )
+   return_value = 0
+
+   WRITE( temp_string, * ) file_id
+   CALL internal_message( 'debug',        &
+                          routine_name // &
+                          ': finalize file definition (file_id=' // TRIM( temp_string ) // ')' )
+
+   !-- Set general no fill, otherwise the performance drops significantly
+   nc_stat = NF90_SET_FILL( file_id, NF90_NOFILL, old_mode )
+
+   IF ( nc_stat == NF90_NOERR )  THEN
+      nc_stat = NF90_ENDDEF( file_id )
+   ENDIF
+
+   IF ( nc_stat /= NF90_NOERR )  THEN
+      return_value = 1
+      CALL internal_message( 'error', routine_name // ': NetCDF error: ' // &
+                                      NF90_STRERROR( nc_stat ) )
+   ENDIF
+#else
+   return_value = 1
+#endif
+
+END SUBROUTINE netcdf4_serial_init_end
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Write variable of different kind into netcdf file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_serial_write_variable(                                 &
+              file_id, var_id, bounds_start, bounds_end, bounds_origin,   &
+              do_output, is_global,                                       &
+              var_int8_0d,   var_int8_1d,   var_int8_2d,   var_int8_3d,   &
+              var_int16_0d,  var_int16_1d,  var_int16_2d,  var_int16_3d,  &
+              var_int32_0d,  var_int32_1d,  var_int32_2d,  var_int32_3d,  &
+              var_intwp_0d,  var_intwp_1d,  var_intwp_2d,  var_intwp_3d,  &
+              var_real32_0d, var_real32_1d, var_real32_2d, var_real32_3d, &
+              var_real64_0d, var_real64_1d, var_real64_2d, var_real64_3d, &
+              var_realwp_0d, var_realwp_1d, var_realwp_2d, var_realwp_3d, &
+              return_value )
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_serial_write_variable'  !< name of this routine
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp)              ::  myid = 0      !< id number of processor element
+   INTEGER(iwp)              ::  nc_stat       !< netcdf return value
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+   INTEGER(iwp), INTENT(IN)  ::  var_id        !< variable ID
+
+   INTEGER(iwp), DIMENSION(:), INTENT(IN)  ::  bounds_origin  !< starting index of each dimension
+   INTEGER(iwp), DIMENSION(:), INTENT(IN)  ::  bounds_end     !< ending index of variable
+   INTEGER(iwp), DIMENSION(:), INTENT(IN)  ::  bounds_start   !< starting index of variable
+   INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  value_count    !< count of values along each dimension to be written
+
+   INTEGER(KIND=1), POINTER, INTENT(IN), OPTIONAL                               ::  var_int8_0d  !< output variable
+   INTEGER(KIND=1), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_int8_1d  !< output variable
+   INTEGER(KIND=1), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_int8_2d  !< output variable
+   INTEGER(KIND=1), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_int8_3d  !< output variable
+
+   INTEGER(KIND=2), POINTER, INTENT(IN), OPTIONAL                               ::  var_int16_0d  !< output variable
+   INTEGER(KIND=2), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_int16_1d  !< output variable
+   INTEGER(KIND=2), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_int16_2d  !< output variable
+   INTEGER(KIND=2), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_int16_3d  !< output variable
+
+   INTEGER(KIND=4), POINTER, INTENT(IN), OPTIONAL                               ::  var_int32_0d  !< output variable
+   INTEGER(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_int32_1d  !< output variable
+   INTEGER(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_int32_2d  !< output variable
+   INTEGER(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_int32_3d  !< output variable
+
+   INTEGER(iwp), POINTER, INTENT(IN), OPTIONAL                               ::  var_intwp_0d  !< output variable
+   INTEGER(iwp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_intwp_1d  !< output variable
+   INTEGER(iwp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_intwp_2d  !< output variable
+   INTEGER(iwp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_intwp_3d  !< output variable
+
+   LOGICAL, INTENT(IN) ::  do_output  !< if false, set count to 0 and do no output
+   LOGICAL, INTENT(IN) ::  is_global  !< true if variable is global (same on all PE)
+
+   REAL(KIND=4), POINTER, INTENT(IN), OPTIONAL                               ::  var_real32_0d  !< output variable
+   REAL(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_real32_1d  !< output variable
+   REAL(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_real32_2d  !< output variable
+   REAL(KIND=4), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_real32_3d  !< output variable
+
+   REAL(KIND=8), POINTER, INTENT(IN), OPTIONAL                               ::  var_real64_0d  !< output variable
+   REAL(KIND=8), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_real64_1d  !< output variable
+   REAL(KIND=8), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_real64_2d  !< output variable
+   REAL(KIND=8), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_real64_3d  !< output variable
+
+   REAL(wp), POINTER, INTENT(IN), OPTIONAL                               ::  var_realwp_0d  !< output variable
+   REAL(wp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:)     ::  var_realwp_1d  !< output variable
+   REAL(wp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:)   ::  var_realwp_2d  !< output variable
+   REAL(wp), POINTER, CONTIGUOUS, INTENT(IN), OPTIONAL, DIMENSION(:,:,:) ::  var_realwp_3d  !< output variable
+
+
+#if defined( __netcdf4 )
+   return_value = 0
+
+#if defined( __parallel )
+   CALL MPI_COMM_RANK( MPI_COMM_WORLD, myid, return_value )
+   IF ( return_value /= 0 )  THEN
+      CALL internal_message( 'error', routine_name // ': MPI_COMM_RANK error' )
+   ENDIF
+#endif
+
+   IF ( myid == 0 )  THEN
+
+      IF ( is_global )  CONTINUE  ! reqired to prevent compiler warning
+
+      WRITE( temp_string, * ) var_id
+      CALL internal_message( 'debug', routine_name // ': write variable ' // TRIM( temp_string ) )
+
+      ALLOCATE( value_count(SIZE( bounds_start )) )
+
+      IF ( do_output ) THEN
+         value_count = bounds_end - bounds_start + 1
+      ELSE
+         value_count = 0
+      END IF
+
+      !-- 8bit integer output
+      IF ( PRESENT( var_int8_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_int8_0d /),       &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int8_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int8_1d,             &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int8_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int8_2d,             &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int8_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int8_3d,             &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      !-- 16bit integer output
+      ELSEIF ( PRESENT( var_int16_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_int16_0d /),      &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int16_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int16_1d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int16_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int16_2d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int16_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int16_3d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      !-- 32bit integer output
+      ELSEIF ( PRESENT( var_int32_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_int32_0d /),      &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int32_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int32_1d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int32_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int32_2d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_int32_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_int32_3d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      !-- working-precision integer output
+      ELSEIF ( PRESENT( var_intwp_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_intwp_0d /),      &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_intwp_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_intwp_1d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_intwp_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_intwp_2d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_intwp_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_intwp_3d,            &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      !-- 32bit real output
+      ELSEIF ( PRESENT( var_real32_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_real32_0d /),     &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_real32_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_real32_1d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_real32_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_real32_2d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_real32_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_real32_3d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      !-- 64bit real output
+      ELSEIF ( PRESENT( var_real64_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_real64_0d /),     &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_real64_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_real64_1d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_real64_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_real64_2d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_real64_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_real64_3d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      !-- working-precision real output
+      ELSEIF ( PRESENT( var_realwp_0d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, (/ var_realwp_0d /),     &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_realwp_1d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_realwp_1d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_realwp_2d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_realwp_2d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSEIF ( PRESENT( var_realwp_3d ) )  THEN
+         nc_stat = NF90_PUT_VAR( file_id, var_id, var_realwp_3d,           &
+                                 start = bounds_start - bounds_origin + 1, &
+                                 count = value_count )
+      ELSE
+         return_value = 1
+         nc_stat = NF90_NOERR
+         WRITE( temp_string, '(": variable_id=",I6,", file_id=",I6)' ) var_id, file_id
+         CALL internal_message( 'error', routine_name // &
+                                TRIM( temp_string ) //   &
+                                ': no output values given' )
+      ENDIF
+
+      !-- Check for errors
+      IF ( nc_stat /= NF90_NOERR )  THEN
+         return_value = 1
+         WRITE( temp_string, * ) 'variable_id=', var_id, '; file_id=', file_id, &
+                                 ', bounds_start=', bounds_start,               &
+                                 ', bounds_end=', bounds_end,                   &
+                                 ', count=', value_count
+         CALL internal_message( 'error', routine_name //                                    &
+                                         ': error while writing ' // TRIM( temp_string ) // &
+                                         ': NetCDF error: ' // NF90_STRERROR( nc_stat ) )
+      ENDIF
+
+   ENDIF
+#else
+   return_value = 1
+#endif
+
+END SUBROUTINE netcdf4_serial_write_variable
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Close netcdf file.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_serial_finalize( file_id, return_value )
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'netcdf4_serial_finalize'  !< name of routine
+
+   INTEGER(iwp), INTENT(IN)  ::  file_id       !< file ID
+   INTEGER(iwp)              ::  nc_stat       !< netcdf return value
+   INTEGER(iwp), INTENT(OUT) ::  return_value  !< return value
+
+
+#if defined( __netcdf4 )
+   WRITE( temp_string, * ) file_id
+   CALL internal_message( 'debug', routine_name // &
+                                   ': close file (file_id=' // TRIM( temp_string ) // ')' )
+
+   nc_stat = NF90_CLOSE( file_id )
+   IF ( nc_stat == NF90_NOERR )  THEN
+      return_value = 0
+   ELSE
+      return_value = 1
+      CALL internal_message( 'error', routine_name // &
+                                      ': NetCDF error: ' // NF90_STRERROR( nc_stat ) )
+   ENDIF
+#else
+   return_value = 1
+#endif
+
+END SUBROUTINE netcdf4_serial_finalize
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Convert data_type string into netcdf data type value.
+!--------------------------------------------------------------------------------------------------!
+FUNCTION get_netcdf_data_type( data_type ) RESULT( return_value )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  data_type  !< requested data type
+
+   CHARACTER(LEN=*), PARAMETER ::  routine_name = 'get_netcdf_data_type'  !< name of this routine
+
+   INTEGER(iwp) ::  return_value  !< netcdf data type
+
+
+   SELECT CASE ( TRIM( data_type ) )
+
+#if defined( __netcdf4 )
+      CASE ( 'char' )
+         return_value = NF90_CHAR
+
+      CASE ( 'int8' )
+         return_value = NF90_BYTE
+
+      CASE ( 'int16' )
+         return_value = NF90_SHORT
+
+      CASE ( 'int32' )
+         return_value = NF90_INT
+
+      CASE ( 'real32' )
+         return_value = NF90_FLOAT
+
+      CASE ( 'real64' )
+         return_value = NF90_DOUBLE
+#endif
+
+      CASE DEFAULT
+         CALL internal_message( 'error', routine_name // &
+                                         ': data type unknown (' // TRIM( data_type ) // ')' )
+         return_value = -1_iwp
+
+   END SELECT
+
+END FUNCTION get_netcdf_data_type
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Message routine writing debug information into the debug file
+!> or creating the error message string.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE internal_message( level, string )
+
+   CHARACTER(LEN=*), INTENT(IN) ::  level   !< message importance level
+   CHARACTER(LEN=*), INTENT(IN) ::  string  !< message string
+
+
+   IF ( TRIM( level ) == 'error' )  THEN
+
+      WRITE( internal_error_message, '(A,A)' ) 'DOM ERROR: ', string
+
+   ELSEIF ( TRIM( level ) == 'debug'  .AND.  print_debug_output )  THEN
+
+      WRITE( debug_output_unit, '(A,A)' ) 'DOM DEBUG: ', string
+      FLUSH( debug_output_unit )
+
+   ENDIF
+
+END SUBROUTINE internal_message
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Return the last created error message.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE netcdf4_serial_get_error_message( error_message )
+
+   CHARACTER(LEN=800), INTENT(OUT) ::  error_message  !< return error message to main program
+
+
+   error_message = internal_error_message
+
+END SUBROUTINE netcdf4_serial_get_error_message
+
+
+END MODULE data_output_netcdf4_serial_module
