Index: palm/trunk/SOURCE/data_output_module.f90
===================================================================
--- palm/trunk/SOURCE/data_output_module.f90	(revision 4116)
+++ palm/trunk/SOURCE/data_output_module.f90	(revision 4123)
@@ -89,5 +89,5 @@
       CHARACTER(LEN=charlen) ::  data_type = ''       !< data type
       CHARACTER(LEN=charlen) ::  name                 !< variable name
-      INTEGER(iwp)           ::  id = 0               !< id within file
+      INTEGER(iwp)           ::  id = -1              !< id within file
       LOGICAL                ::  is_global = .FALSE.  !< true if global variable
       CHARACTER(LEN=charlen), DIMENSION(:), ALLOCATABLE ::  dimension_names  !< list of dimension names
@@ -99,8 +99,8 @@
       CHARACTER(LEN=charlen) ::  data_type = ''       !< data type
       CHARACTER(LEN=charlen) ::  name                 !< dimension name
-      INTEGER(iwp)           ::  id = 0               !< dimension id within file
+      INTEGER(iwp)           ::  id = -1              !< 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
+      INTEGER(iwp)           ::  var_id = -1          !< associated variable id within file
       LOGICAL                ::  is_masked = .FALSE.  !< true if masked
       INTEGER(iwp),    DIMENSION(2)              ::  bounds                !< lower and upper bound of dimension
@@ -127,5 +127,5 @@
       CHARACTER(LEN=charlen) ::  format = ''        !< file format
       CHARACTER(LEN=charlen) ::  name = ''          !< file name
-      INTEGER(iwp)           ::  id = 0             !< id of file
+      INTEGER(iwp)           ::  id = -1            !< id of file
       LOGICAL                ::  is_init = .FALSE.  !< true if initialized
       TYPE(attribute_type), DIMENSION(:), ALLOCATABLE ::  attributes  !< list of attributes
@@ -965,5 +965,5 @@
                ALLOCATE( dimension%masked_values_realwp(0:dimension%length_mask-1) )
                j = 0
-               DO  i = 0, dimension%length-1
+               DO  i = dimension%bounds(1), dimension%bounds(2)   !> @todo change loop also for other data types
                   IF ( dimension%mask(i) )  THEN
                      dimension%masked_values_realwp(j) = dimension%values_realwp(i)
@@ -1085,4 +1085,21 @@
 ! ------------
 !> Add variable to database.
+!> Example call:
+!>   dom_def_var( filename =  'DATA_OUTPUT_3D', &
+!>                name = 'u', &
+!>                dimension_names = (/'x   ', 'y   ', 'z   ', 'time'/), &
+!>                output_type = 'real32' )
+!> @note The order of dimensions must match in reversed order to the dimensions of the
+!>       corresponding variable array. The last given dimension can also be non-existent within the
+!>       variable array if at any given call of 'dom_write_var' for this variable, the last
+!>       dimension has only a single index.
+!>       Hence, the array 'u' must be allocated with dimension 'x' as its last dimension, preceded
+!>       by 'y', then 'z', and 'time' being the first dimension. If at any given write statement,
+!>       only a single index of dimension 'time' is to be written, the dimension can be non-present
+!>       in the variable array leaving dimension 'z' as the first dimension.
+!>       So, the variable array needs to be allocated like either:
+!>          ALLOCATE( u(<time>,<z>,<y>,<x>) )
+!>       or
+!>          ALLOCATE( u(<z>,<y>,<x>) )
 !--------------------------------------------------------------------------------------------------!
 FUNCTION dom_def_var( filename, name, dimension_names, output_type, is_global ) &
@@ -1123,4 +1140,5 @@
 
    variable%dimension_names = dimension_names
+   variable%dimension_ids = -1
    variable%data_type = TRIM( output_type )
 
@@ -2578,6 +2596,12 @@
 !>                  name = 'u', &
 !>                  var_real64_3d = u, &
-!>                  bounds_start = (/nzb, nys, nxl/), &
-!>                  bounds_end = (/nzt, nyn, nxr/)  )
+!>                  bounds_start = (/nxl, nys, nzb, time_step/), &
+!>                  bounds_end = (/nxr, nyn, nzt, time_step/)  )
+!> @note The order of dimension bounds must match to the order of dimensions given in call
+!>       'dom_def_var'. I.e., the corresponding variable definition should be like:
+!>          dom_def_var( filename =  'DATA_OUTPUT_3D', &
+!>                       name = 'u', &
+!>                       dimension_names = (/'x   ', 'y   ', 'z   ', 'time'/), &
+!>                       output_type = <desired-output-type> )
 !--------------------------------------------------------------------------------------------------!
 FUNCTION dom_write_var( filename, name, bounds_start, bounds_end,       &
@@ -2606,10 +2630,10 @@
    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 msked var at each dim
-   INTEGER(iwp), DIMENSION(:),   ALLOCATABLE ::  bounds_start_new  !< start index (lower bound) of msked var at each dim
-   INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  masked_indices    !< dummy list holding all masked indices along a dimension
+   INTEGER(iwp), DIMENSION(:),   INTENT(IN)  ::  bounds_end             !< end index per dimension of variable
+   INTEGER(iwp), DIMENSION(:),   INTENT(IN)  ::  bounds_start           !< start index per dimension of variable
+   INTEGER(iwp), DIMENSION(:),   ALLOCATABLE ::  bounds_origin          !< first index of each dimension
+   INTEGER(iwp), DIMENSION(:),   ALLOCATABLE ::  bounds_start_internal  !< start index per dim. for output after masking
+   INTEGER(iwp), DIMENSION(:),   ALLOCATABLE ::  value_counts           !< count of indices to be written per dimension
+   INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  masked_indices         !< list containing all output indices along a dimension
 
    LOGICAL ::  do_output  !< true if any data lies within given range of masked dimension
@@ -2745,11 +2769,7 @@
 
       !-- 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 ))   )
+      ALLOCATE( bounds_origin(SIZE( dimension_list )) )
+      ALLOCATE( bounds_start_internal(SIZE( dimension_list )) )
+      ALLOCATE( value_counts(SIZE( dimension_list )) )
 
       WRITE( temp_string, * ) bounds_start
@@ -2765,18 +2785,20 @@
 
       !-- 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 get_masked_indices_and_masked_dimension_bounds( dimension_list,                  &
+              bounds_start, bounds_end, bounds_start_internal, value_counts, bounds_origin, &
+              masked_indices )
+
+      do_output = .NOT. ANY( value_counts == 0 )
+
+      WRITE( temp_string, * ) bounds_start_internal
       CALL internal_message( 'debug', routine_name //                     &
                                       ': file "' // TRIM( filename ) //   &
                                        '": variable "' // TRIM( name ) // &
-                                       '": bounds_start_new =' // TRIM( temp_string ) )
-      WRITE( temp_string, * ) bounds_end_new
+                                       '": bounds_start_internal =' // TRIM( temp_string ) )
+      WRITE( temp_string, * ) value_counts
       CALL internal_message( 'debug', routine_name //                     &
                                       ': file "' // TRIM( filename ) //   &
                                        '": variable "' // TRIM( name ) // &
-                                       '": bounds_end_new =' // TRIM( temp_string ) )
+                                       '": value_counts =' // TRIM( temp_string ) )
 
       !-- Mask and resort variable
@@ -2786,8 +2808,8 @@
       ELSEIF ( PRESENT( var_int8_1d ) )  THEN
          IF ( do_output ) THEN
-            ALLOCATE( var_int8_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            ALLOCATE( var_int8_1d_resorted(0:value_counts(1)-1) )
             !$OMP PARALLEL PRIVATE (i)
             !$OMP DO
-            DO  i = bounds_start_new(1), bounds_end_new(1)
+            DO  i = 0, value_counts(1) - 1
                var_int8_1d_resorted(i) = var_int8_1d(masked_indices(1,i))
             ENDDO
@@ -2800,10 +2822,10 @@
       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)) )
+            ALLOCATE( var_int8_2d_resorted(0:value_counts(1)-1, &
+                                           0:value_counts(2)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
                   var_int8_2d_resorted(i,j) = var_int8_2d(masked_indices(2,j), &
                                                           masked_indices(1,i)  )
@@ -2818,12 +2840,12 @@
       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)) )
+            ALLOCATE( var_int8_3d_resorted(0:value_counts(1)-1, &
+                                           0:value_counts(2)-1, &
+                                           0:value_counts(3)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
+                  DO  k = 0, value_counts(3) - 1
                      var_int8_3d_resorted(i,j,k) = var_int8_3d(masked_indices(3,k), &
                                                                masked_indices(2,j), &
@@ -2844,8 +2866,8 @@
       ELSEIF ( PRESENT( var_int16_1d ) )  THEN
          IF ( do_output ) THEN
-            ALLOCATE( var_int16_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            ALLOCATE( var_int16_1d_resorted(0:value_counts(1)-1) )
             !$OMP PARALLEL PRIVATE (i)
             !$OMP DO
-            DO  i = bounds_start_new(1), bounds_end_new(1)
+            DO  i = 0, value_counts(1) - 1
                var_int16_1d_resorted(i) = var_int16_1d(masked_indices(1,i))
             ENDDO
@@ -2858,10 +2880,10 @@
       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)) )
+            ALLOCATE( var_int16_2d_resorted(0:value_counts(1)-1, &
+                                            0:value_counts(2)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
                   var_int16_2d_resorted(i,j) = var_int16_2d(masked_indices(2,j), &
                                                             masked_indices(1,i))
@@ -2876,12 +2898,12 @@
       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)) )
+            ALLOCATE( var_int16_3d_resorted(0:value_counts(1)-1, &
+                                            0:value_counts(2)-1, &
+                                            0:value_counts(3)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
+                  DO  k = 0, value_counts(3) - 1
                      var_int16_3d_resorted(i,j,k) = var_int16_3d(masked_indices(3,k), &
                                                                  masked_indices(2,j), &
@@ -2902,8 +2924,8 @@
       ELSEIF ( PRESENT( var_int32_1d ) )  THEN
          IF ( do_output ) THEN
-            ALLOCATE( var_int32_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            ALLOCATE( var_int32_1d_resorted(0:value_counts(1)-1) )
             !$OMP PARALLEL PRIVATE (i)
             !$OMP DO
-            DO  i = bounds_start_new(1), bounds_end_new(1)
+            DO  i = 0, value_counts(1) - 1
                var_int32_1d_resorted(i) = var_int32_1d(masked_indices(1,i))
             ENDDO
@@ -2916,10 +2938,10 @@
       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)) )
+            ALLOCATE( var_int32_2d_resorted(0:value_counts(1)-1, &
+                                            0:value_counts(2)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
                   var_int32_2d_resorted(i,j) = var_int32_2d(masked_indices(2,j), &
                                                             masked_indices(1,i)  )
@@ -2934,12 +2956,12 @@
       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)) )
+            ALLOCATE( var_int32_3d_resorted(0:value_counts(1)-1, &
+                                            0:value_counts(2)-1, &
+                                            0:value_counts(3)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
+                  DO  k = 0, value_counts(3) - 1
                      var_int32_3d_resorted(i,j,k) = var_int32_3d(masked_indices(3,k), &
                                                                  masked_indices(2,j), &
@@ -2960,8 +2982,8 @@
       ELSEIF ( PRESENT( var_intwp_1d ) )  THEN
          IF ( do_output ) THEN
-            ALLOCATE( var_intwp_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            ALLOCATE( var_intwp_1d_resorted(0:value_counts(1)-1) )
             !$OMP PARALLEL PRIVATE (i)
             !$OMP DO
-            DO  i = bounds_start_new(1), bounds_end_new(1)
+            DO  i = 0, value_counts(1) - 1
                var_intwp_1d_resorted(i) = var_intwp_1d(masked_indices(1,i))
             ENDDO
@@ -2974,10 +2996,10 @@
       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)) )
+            ALLOCATE( var_intwp_2d_resorted(0:value_counts(1)-1, &
+                                            0:value_counts(2)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
                   var_intwp_2d_resorted(i,j) = var_intwp_2d(masked_indices(2,j), &
                                                             masked_indices(1,i)  )
@@ -2992,12 +3014,12 @@
       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)) )
+            ALLOCATE( var_intwp_3d_resorted(0:value_counts(1)-1, &
+                                            0:value_counts(2)-1, &
+                                            0:value_counts(3)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
+                  DO  k = 0, value_counts(3) - 1
                      var_intwp_3d_resorted(i,j,k) = var_intwp_3d(masked_indices(3,k), &
                                                                  masked_indices(2,j), &
@@ -3018,8 +3040,8 @@
       ELSEIF ( PRESENT( var_real32_1d ) )  THEN
          IF ( do_output ) THEN
-            ALLOCATE( var_real32_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            ALLOCATE( var_real32_1d_resorted(0:value_counts(1)-1) )
             !$OMP PARALLEL PRIVATE (i)
             !$OMP DO
-            DO  i = bounds_start_new(1), bounds_end_new(1)
+            DO  i = 0, value_counts(1) - 1
                var_real32_1d_resorted(i) = var_real32_1d(masked_indices(1,i))
             ENDDO
@@ -3032,10 +3054,10 @@
       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)) )
+            ALLOCATE( var_real32_2d_resorted(0:value_counts(1)-1, &
+                                             0:value_counts(2)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
                   var_real32_2d_resorted(i,j) = var_real32_2d(masked_indices(2,j), &
                                                               masked_indices(1,i)  )
@@ -3050,12 +3072,12 @@
       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)) )
+            ALLOCATE( var_real32_3d_resorted(0:value_counts(1)-1, &
+                                             0:value_counts(2)-1, &
+                                             0:value_counts(3)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
+                  DO  k = 0, value_counts(3) - 1
                      var_real32_3d_resorted(i,j,k) = var_real32_3d(masked_indices(3,k), &
                                                                    masked_indices(2,j), &
@@ -3076,8 +3098,8 @@
       ELSEIF ( PRESENT( var_real64_1d ) )  THEN
          IF ( do_output ) THEN
-            ALLOCATE( var_real64_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            ALLOCATE( var_real64_1d_resorted(0:value_counts(1)-1) )
             !$OMP PARALLEL PRIVATE (i)
             !$OMP DO
-            DO  i = bounds_start_new(1), bounds_end_new(1)
+            DO  i = 0, value_counts(1) - 1
                var_real64_1d_resorted(i) = var_real64_1d(masked_indices(1,i))
             ENDDO
@@ -3090,10 +3112,10 @@
       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)) )
+            ALLOCATE( var_real64_2d_resorted(0:value_counts(1)-1, &
+                                             0:value_counts(2)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
                   var_real64_2d_resorted(i,j) = var_real64_2d(masked_indices(2,j), &
                                                               masked_indices(1,i)  )
@@ -3108,12 +3130,12 @@
       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)) )
+            ALLOCATE( var_real64_3d_resorted(0:value_counts(1)-1, &
+                                             0:value_counts(2)-1, &
+                                             0:value_counts(3)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
+                  DO  k = 0, value_counts(3) - 1
                      var_real64_3d_resorted(i,j,k) = var_real64_3d(masked_indices(3,k), &
                                                                    masked_indices(2,j), &
@@ -3134,8 +3156,8 @@
       ELSEIF ( PRESENT( var_realwp_1d ) )  THEN
          IF ( do_output ) THEN
-            ALLOCATE( var_realwp_1d_resorted(bounds_start_new(1):bounds_end_new(1)) )
+            ALLOCATE( var_realwp_1d_resorted(0:value_counts(1)-1) )
             !$OMP PARALLEL PRIVATE (i)
             !$OMP DO
-            DO  i = bounds_start_new(1), bounds_end_new(1)
+            DO  i = 0, value_counts(1) - 1
                var_realwp_1d_resorted(i) = var_realwp_1d(masked_indices(1,i))
             ENDDO
@@ -3148,10 +3170,10 @@
       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)) )
+            ALLOCATE( var_realwp_2d_resorted(0:value_counts(1)-1, &
+                                             0:value_counts(2)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
                   var_realwp_2d_resorted(i,j) = var_realwp_2d(masked_indices(2,j), &
                                                               masked_indices(1,i)  )
@@ -3166,12 +3188,12 @@
       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)) )
+            ALLOCATE( var_realwp_3d_resorted(0:value_counts(1)-1, &
+                                             0:value_counts(2)-1, &
+                                             0:value_counts(3)-1) )
             !$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)
+            DO  i = 0, value_counts(1) - 1
+               DO  j = 0, value_counts(2) - 1
+                  DO  k = 0, value_counts(3) - 1
                      var_realwp_3d_resorted(i,j,k) = var_realwp_3d(masked_indices(3,k), &
                                                                    masked_indices(2,j), &
@@ -3207,120 +3229,120 @@
             !-- 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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int8_0d=var_int8_0d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int8_1d=var_int8_1d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int8_2d=var_int8_2d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int8_3d=var_int8_3d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int16_0d=var_int16_0d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int16_1d=var_int16_1d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int16_2d=var_int16_2d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int16_3d=var_int16_3d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int32_0d=var_int32_0d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int32_1d=var_int32_1d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int32_2d=var_int32_2d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int32_3d=var_int32_3d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_intwp_0d=var_intwp_0d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_intwp_1d=var_intwp_1d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_intwp_2d=var_intwp_2d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_intwp_3d=var_intwp_3d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real32_0d=var_real32_0d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real32_1d=var_real32_1d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real32_2d=var_real32_2d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real32_3d=var_real32_3d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real64_0d=var_real64_0d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real64_1d=var_real64_1d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real64_2d=var_real64_2d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real64_3d=var_real64_3d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_realwp_0d=var_realwp_0d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_realwp_1d=var_realwp_1d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_realwp_2d=var_realwp_2d_pointer, return_value=output_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, &
+               CALL binary_write_variable( file_id, var_id,                           &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_realwp_3d=var_realwp_3d_pointer, return_value=output_return_value )
             ELSE
@@ -3336,120 +3358,120 @@
             !-- 8bit integer output
             IF ( PRESENT( var_int8_0d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int8_0d=var_int8_0d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_int8_1d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int8_1d=var_int8_1d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_int8_2d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int8_2d=var_int8_2d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_int8_3d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int8_3d=var_int8_3d_pointer, return_value=output_return_value )
             !-- 16bit integer output
             ELSEIF ( PRESENT( var_int16_0d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int16_0d=var_int16_0d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_int16_1d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int16_1d=var_int16_1d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_int16_2d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int16_2d=var_int16_2d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_int16_3d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int16_3d=var_int16_3d_pointer, return_value=output_return_value )
             !-- 32bit integer output
             ELSEIF ( PRESENT( var_int32_0d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int32_0d=var_int32_0d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_int32_1d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int32_1d=var_int32_1d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_int32_2d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int32_2d=var_int32_2d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_int32_3d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_int32_3d=var_int32_3d_pointer, return_value=output_return_value )
             !-- working-precision integer output
             ELSEIF ( PRESENT( var_intwp_0d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_intwp_0d=var_intwp_0d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_intwp_1d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_intwp_1d=var_intwp_1d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_intwp_2d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_intwp_2d=var_intwp_2d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_intwp_3d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_intwp_3d=var_intwp_3d_pointer, return_value=output_return_value )
             !-- 32bit real output
             ELSEIF ( PRESENT( var_real32_0d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real32_0d=var_real32_0d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_real32_1d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real32_1d=var_real32_1d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_real32_2d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real32_2d=var_real32_2d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_real32_3d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real32_3d=var_real32_3d_pointer, return_value=output_return_value )
             !-- 64bit real output
             ELSEIF ( PRESENT( var_real64_0d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real64_0d=var_real64_0d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_real64_1d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real64_1d=var_real64_1d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_real64_2d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real64_2d=var_real64_2d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_real64_3d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_real64_3d=var_real64_3d_pointer, return_value=output_return_value )
             !-- working-precision real output
             ELSEIF ( PRESENT( var_realwp_0d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_realwp_0d=var_realwp_0d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_realwp_1d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_realwp_1d=var_realwp_1d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_realwp_2d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_realwp_2d=var_realwp_2d_pointer, return_value=output_return_value )
             ELSEIF ( PRESENT( var_realwp_3d ) )  THEN
-               CALL netcdf4_write_variable( file_id, var_id,                              &
-                       bounds_start_new, bounds_end_new, bounds_dim_start, do_output, is_global, &
+               CALL netcdf4_write_variable( file_id, var_id,                          &
+                       bounds_start_internal, value_counts, bounds_origin, is_global, &
                        var_realwp_3d=var_realwp_3d_pointer, return_value=output_return_value )
             ELSE
@@ -3608,11 +3630,13 @@
 ! ------------
 !> 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 )
+!> 'bounds_end'). Return the masked indices ('masked_indices') of the dimensions, the first index
+!> of the masked dimensions containing these indices ('bounds_masked_start'), the count of masked
+!> indices within given bounds ('value_counts') and the origin index of each dimension
+!> ('bounds_origin'). If, for any dimension, no masked index lies within the given bounds, counts,
+!> starts and origins are set to zero for all dimensions.
+!--------------------------------------------------------------------------------------------------!
+SUBROUTINE get_masked_indices_and_masked_dimension_bounds(                             &
+              dimensions, bounds_start, bounds_end, bounds_masked_start, value_counts, &
+              bounds_origin, masked_indices )
 
    ! CHARACTER(LEN=*), PARAMETER ::  routine_name = 'get_masked_indices_and_masked_dimension_bounds'  !< name of routine
@@ -3623,17 +3647,15 @@
 
    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(OUT) ::  bounds_origin        !< first index of each dimension, 0 if dimension is masked
    INTEGER(iwp), DIMENSION(:), INTENT(IN)  ::  bounds_start         !< lower bounds to be searched in
+   INTEGER(iwp), DIMENSION(:), INTENT(OUT) ::  value_counts         !< count of indices per dimension to be output
 
    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 )) )
+   ALLOCATE( masked_indices(SIZE( dimensions ),0:MAXVAL( bounds_end - bounds_start + 1 )) )
    masked_indices = -HUGE( 0_iwp )
 
@@ -3643,23 +3665,23 @@
       IF ( dimensions(d)%is_masked )  THEN
 
-         bounds_masked_end(d) = HUGE( 0_iwp )
+         bounds_origin(d) = 0
+
          bounds_masked_start(d) = -HUGE( 0_iwp )
 
          !-- Find number of masked values within given variable bounds
-         j = 0
+         value_counts(d) = 0
          DO  i = LBOUND( dimensions(d)%masked_indices, DIM=1 ), &
                  UBOUND( dimensions(d)%masked_indices, DIM=1 )
 
-            !-- Is masked position within given bounds?
+            !-- Is masked index 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 masked index
+               masked_indices(d,value_counts(d)) = dimensions(d)%masked_indices(i)
+               value_counts(d) = value_counts(d) + 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
@@ -3667,9 +3689,9 @@
          ENDDO
 
-         !-- Set masked bounds to zero if no masked data lies within bounds
-         IF ( j == 0 )  THEN
+         !-- Set masked bounds to zero if no masked index lies within bounds
+         IF ( value_counts(d) == 0 )  THEN
+            bounds_origin(:) = 0
             bounds_masked_start(:) = 0_iwp
-            bounds_masked_end(:) = 0_iwp
-            mask_in_bounds = .FALSE.
+            value_counts(:) = 0_iwp
             EXIT
          ENDIF
@@ -3677,9 +3699,11 @@
       ELSE
 
-         !-- If dimension is not masked, tread all data within bounds as masked
+         !-- If dimension is not masked, save all indices within bounds for output
+         bounds_origin(d) = dimensions(d)%bounds(1)
          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
+         value_counts(d) = bounds_end(d) - bounds_start(d) + 1
+
+         DO  i = 0, value_counts(d) - 1
+            masked_indices(d,i) = bounds_start(d) + i
          ENDDO
 
