Index: palm/trunk/SOURCE/netcdf.f90
===================================================================
--- palm/trunk/SOURCE/netcdf.f90	(revision 392)
+++ palm/trunk/SOURCE/netcdf.f90	(revision 410)
@@ -8,4 +8,8 @@
 ! ------------------
 ! 
+!
+! Branch revisions:
+! -----------------
+! masked data output
 !
 ! Former revisions:
@@ -94,9 +98,10 @@
     CHARACTER (LEN=2000)           ::  var_list, var_list_old
 
-    INTEGER ::  av, i, id_x, id_y, id_z, j, ns, ns_old, nz_old
+    INTEGER ::  av, file_id, i, id_x, id_y, id_z, j, ns, ns_old, nz_old
 
     INTEGER, DIMENSION(1) ::  id_dim_time_old, id_dim_x_yz_old,  &
                               id_dim_y_xz_old, id_dim_zu_sp_old, &
-                              id_dim_zu_xy_old, id_dim_zu_3d_old
+                              id_dim_zu_xy_old, id_dim_zu_3d_old, &
+                              id_dim_zu_mask_old
 
     LOGICAL ::  found
@@ -108,5 +113,7 @@
     REAL, DIMENSION(1) ::  last_time_coordinate
 
-    REAL, DIMENSION(:), ALLOCATABLE ::  netcdf_data
+    REAL, DIMENSION(:), ALLOCATABLE   ::  netcdf_data
+    REAL, DIMENSION(:,:), ALLOCATABLE ::  netcdf_data_2d
+
 
 
@@ -161,4 +168,86 @@
              CASE ( 'prt' )
                 nc_precision(8) = j
+             CASE ( 'masks' )
+                nc_precision(11:50) = j
+             CASE ( 'mask01' )
+                nc_precision(11) = j
+             CASE ( 'mask02' )
+                nc_precision(12) = j
+             CASE ( 'mask03' )
+                nc_precision(13) = j
+             CASE ( 'mask04' )
+                nc_precision(14) = j
+             CASE ( 'mask05' )
+                nc_precision(15) = j
+             CASE ( 'mask06' )
+                nc_precision(16) = j
+             CASE ( 'mask07' )
+                nc_precision(17) = j
+             CASE ( 'mask08' )
+                nc_precision(18) = j
+             CASE ( 'mask09' )
+                nc_precision(19) = j
+             CASE ( 'mask10' )
+                nc_precision(20) = j
+             CASE ( 'mask11' )
+                nc_precision(21) = j
+             CASE ( 'mask12' )
+                nc_precision(22) = j
+             CASE ( 'mask13' )
+                nc_precision(23) = j
+             CASE ( 'mask14' )
+                nc_precision(24) = j
+             CASE ( 'mask15' )
+                nc_precision(25) = j
+             CASE ( 'mask16' )
+                nc_precision(26) = j
+             CASE ( 'mask17' )
+                nc_precision(27) = j
+             CASE ( 'mask18' )
+                nc_precision(28) = j
+             CASE ( 'mask19' )
+                nc_precision(29) = j
+             CASE ( 'mask20' )
+                nc_precision(30) = j
+             CASE ( 'maskav01' )
+                nc_precision(31) = j
+             CASE ( 'maskav02' )
+                nc_precision(32) = j
+             CASE ( 'maskav03' )
+                nc_precision(33) = j
+             CASE ( 'maskav04' )
+                nc_precision(34) = j
+             CASE ( 'maskav05' )
+                nc_precision(35) = j
+             CASE ( 'maskav06' )
+                nc_precision(36) = j
+             CASE ( 'maskav07' )
+                nc_precision(37) = j
+             CASE ( 'maskav08' )
+                nc_precision(38) = j
+             CASE ( 'maskav09' )
+                nc_precision(39) = j
+             CASE ( 'maskav10' )
+                nc_precision(40) = j
+             CASE ( 'maskav11' )
+                nc_precision(41) = j
+             CASE ( 'maskav12' )
+                nc_precision(42) = j
+             CASE ( 'maskav13' )
+                nc_precision(43) = j
+             CASE ( 'maskav14' )
+                nc_precision(44) = j
+             CASE ( 'maskav15' )
+                nc_precision(45) = j
+             CASE ( 'maskav16' )
+                nc_precision(46) = j
+             CASE ( 'maskav17' )
+                nc_precision(47) = j
+             CASE ( 'maskav18' )
+                nc_precision(48) = j
+             CASE ( 'maskav19' )
+                nc_precision(49) = j
+             CASE ( 'maskav20' )
+                nc_precision(50) = j
              CASE ( 'all' )
                 nc_precision    = j
@@ -173,5 +262,5 @@
 
           i = i + 1
-          IF ( i > 10 )  EXIT
+          IF ( i > 50 )  EXIT
        ENDDO
 
@@ -191,7 +280,563 @@
 
 !
-!-- Select the mode to be processed. Possibilities are 3d, xy, xz, yz,
+!-- Select the mode to be processed. Possibilities are 3d, mask, xy, xz, yz, 
 !-- pr and ts.
     SELECT CASE ( mode )
+
+       CASE ( 'ma_new' )
+
+!
+!--       decompose actual parameter file_id (=formal parameter av) into 
+!--       mid and av
+          file_id = av
+          IF ( file_id <= 140 )  THEN
+             mid = file_id - 120
+             av = 0
+          ELSE
+             mid = file_id - 140
+             av = 1
+          ENDIF
+
+!
+!--       Define some global attributes of the dataset
+          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
+               'Conventions', 'COARDS' )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          IF ( av == 0 )  THEN
+             time_average_text = ' '
+          ELSE
+             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
+                                                            averaging_interval
+          ENDIF
+          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'title', &
+                                  TRIM( run_description_header ) //    &
+                                  TRIM( time_average_text ) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+          IF ( av == 1 )  THEN
+             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
+             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
+                                     'time_avg', TRIM( time_average_text ) )
+             CALL handle_netcdf_error( 'netcdf', 9998 )
+          ENDIF
+
+!
+!--       Define time coordinate for volume data (unlimited dimension)
+          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'time', NF90_UNLIMITED, &
+                                  id_dim_time_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'time', NF90_DOUBLE, &
+                                  id_dim_time_mask(mid,av), &
+                                  id_var_time_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
+                                  id_var_time_mask(mid,av), 'units', &
+                                  'seconds')
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+!
+!--       Define spatial dimensions and coordinates:
+!--       Define vertical coordinate grid (zu grid)
+          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'zu_3d', &
+                                  mask_size(mid,3), id_dim_zu_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zu_3d', NF90_DOUBLE, &
+                                  id_dim_zu_mask(mid,av), &
+                                  id_var_zu_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
+                                  'units', 'meters' )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+!
+!--       Define vertical coordinate grid (zw grid)
+          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'zw_3d', &
+                                  mask_size(mid,3), id_dim_zw_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zw_3d', NF90_DOUBLE, &
+                                  id_dim_zw_mask(mid,av), &
+                                  id_var_zw_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
+                                  'units', 'meters' )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+!
+!--       Define x-axis (for scalar position)
+          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'x', &
+                                  mask_size(mid,1), id_dim_x_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'x', NF90_DOUBLE, &
+                                  id_dim_x_mask(mid,av), id_var_x_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_x_mask(mid,av), &
+                                  'units', 'meters' )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+!
+!--       Define x-axis (for u position)
+          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'xu', &
+                                  mask_size(mid,1), id_dim_xu_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'xu', NF90_DOUBLE, &
+                                  id_dim_xu_mask(mid,av), &
+                                  id_var_xu_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_xu_mask(mid,av), &
+                                  'units', 'meters' )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+!
+!--       Define y-axis (for scalar position)
+          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'y', &
+                                  mask_size(mid,2), id_dim_y_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'y', NF90_DOUBLE, &
+                                  id_dim_y_mask(mid,av), id_var_y_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_y_mask(mid,av), &
+                                  'units', 'meters' )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+!
+!--       Define y-axis (for v position)
+          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'yv', &
+                                  mask_size(mid,2), id_dim_yv_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'yv', NF90_DOUBLE, &
+                                  id_dim_yv_mask(mid,av), &
+                                  id_var_yv_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_yv_mask(mid,av), &
+                                  'units', 'meters' )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+!
+!--       In case of non-flat topography define 2d-arrays containing the height
+!--       informations
+          IF ( TRIM( topography ) /= 'flat' )  THEN
+!
+!--          Define zusi = zu(nzb_s_inner)
+             nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zusi', NF90_DOUBLE, &
+                                     (/ id_dim_x_mask(mid,av),    &
+                                        id_dim_y_mask(mid,av) /), &
+                                     id_var_zusi_mask(mid,av) )
+             CALL handle_netcdf_error( 'netcdf', 9998 )
+             
+             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
+                                     id_var_zusi_mask(mid,av), &
+                                     'units', 'meters' )
+             CALL handle_netcdf_error( 'netcdf', 9998 )
+             
+             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
+                                     id_var_zusi_mask(mid,av), &
+                                     'long_name', 'zu(nzb_s_inner)' )
+             CALL handle_netcdf_error( 'netcdf', 9998 )
+
+!             
+!--          Define zwwi = zw(nzb_w_inner)
+             nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zwwi', NF90_DOUBLE, &
+                                     (/ id_dim_x_mask(mid,av),    &
+                                        id_dim_y_mask(mid,av) /), &
+                                     id_var_zwwi_mask(mid,av) )
+             CALL handle_netcdf_error( 'netcdf', 9998 )
+             
+             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
+                                     id_var_zwwi_mask(mid,av), &
+                                     'units', 'meters' )
+             CALL handle_netcdf_error( 'netcdf', 9998 )
+             
+             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
+                                     id_var_zwwi_mask(mid,av), &
+                                     'long_name', 'zw(nzb_w_inner)' )
+             CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          ENDIF             
+
+
+!
+!--       Define the variables
+          var_list = ';'
+          i = 1
+
+          DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
+
+!
+!--          Check for the grid
+             found = .TRUE.
+             SELECT CASE ( domask(mid,av,i) )
+!
+!--             Most variables are defined on the scalar grid
+                CASE ( 'e', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c', 'ql_v', &
+                       'ql_vp', 'qv', 'rho', 's', 'sa', 'vpt' )
+
+                   grid_x = 'x'
+                   grid_y = 'y'
+                   grid_z = 'zu'
+!
+!--             u grid
+                CASE ( 'u' )
+
+                   grid_x = 'xu'
+                   grid_y = 'y'
+                   grid_z = 'zu'
+!
+!--             v grid
+                CASE ( 'v' )
+
+                   grid_x = 'x'
+                   grid_y = 'yv'
+                   grid_z = 'zu'
+!
+!--             w grid
+                CASE ( 'w' )
+
+                   grid_x = 'x'
+                   grid_y = 'y'
+                   grid_z = 'zw'
+
+                CASE DEFAULT
+!
+!--                Check for user-defined quantities
+                   CALL user_define_netcdf_grid( domask(mid,av,i), found, &
+                                                 grid_x, grid_y, grid_z )
+
+             END SELECT
+
+!
+!--          Select the respective dimension ids
+             IF ( grid_x == 'x' )  THEN
+                id_x = id_dim_x_mask(mid,av)
+             ELSEIF ( grid_x == 'xu' )  THEN
+                id_x = id_dim_xu_mask(mid,av)
+             ENDIF
+
+             IF ( grid_y == 'y' )  THEN
+                id_y = id_dim_y_mask(mid,av)
+             ELSEIF ( grid_y == 'yv' )  THEN
+                id_y = id_dim_yv_mask(mid,av)
+             ENDIF
+
+             IF ( grid_z == 'zu' )  THEN
+                id_z = id_dim_zu_mask(mid,av)
+             ELSEIF ( grid_z == 'zw' )  THEN
+                id_z = id_dim_zw_mask(mid,av)
+             ENDIF
+
+!
+!--          Define the grid
+             nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), domask(mid,av,i), &
+                                     nc_precision(10+mid+av*20),            &
+                                     (/ id_x, id_y, id_z,                   &
+                                        id_dim_time_mask(mid,av) /),        &
+                                     id_var_domask(mid,av,i) )
+
+             IF ( .NOT. found )  THEN
+                WRITE ( message_string, * ) 'no grid defined for', &
+                     ' variable ', TRIM( domask(mid,av,i) )
+                CALL message( 'define_netcdf_header', 'PA9998', 0, 1, 0, 6, 0 )
+             ENDIF
+
+             var_list = TRIM( var_list ) // TRIM( domask(mid,av,i) ) // ';'
+
+             CALL handle_netcdf_error( 'netcdf', 9998 )
+!
+!--          Store the 'real' name of the variable (with *, for example)
+!--          in the long_name attribute. This is evaluated by Ferret,
+!--          for example.
+             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av),     &
+                                     id_var_domask(mid,av,i), &
+                                     'long_name', domask(mid,av,i) )
+             CALL handle_netcdf_error( 'netcdf', 9998 )
+!
+!--          Define the variable's unit
+             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
+                                     id_var_domask(mid,av,i), &
+                                     'units', TRIM( domask_unit(mid,av,i) ) )
+             CALL handle_netcdf_error( 'netcdf', 9998 )
+
+             i = i + 1
+
+          ENDDO
+
+!
+!--       No arrays to output
+          IF ( i == 1 )  RETURN
+
+!
+!--       Write the list of variables as global attribute (this is used by
+!--       restart runs and by combine_plot_fields)
+          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
+                                  'VAR_LIST', var_list )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+!
+!--       Leave NetCDF define mode
+          nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+!
+!--       Write data for x (shifted by +dx/2) and xu axis
+          ALLOCATE( netcdf_data(mask_size(mid,1)) )
+
+          netcdf_data = ( mask_i_global(mid,:mask_size(mid,1)) + 0.5 ) * dx
+
+          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_x_mask(mid,av), &
+                                  netcdf_data, start = (/ 1 /),               &
+                                  count = (/ mask_size(mid,1) /) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          netcdf_data = mask_i_global(mid,:mask_size(mid,1)) * dx
+
+          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_xu_mask(mid,av),&
+                                  netcdf_data, start = (/ 1 /),               &
+                                  count = (/ mask_size(mid,1) /) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          DEALLOCATE( netcdf_data )
+
+!
+!--       Write data for y (shifted by +dy/2) and yv axis
+          ALLOCATE( netcdf_data(mask_size(mid,2)) )
+
+          netcdf_data = ( mask_j_global(mid,:mask_size(mid,2)) + 0.5 ) * dy
+
+          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_y_mask(mid,av), &
+                                  netcdf_data, start = (/ 1 /),               &
+                                  count = (/ mask_size(mid,2) /))
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          netcdf_data = mask_j_global(mid,:mask_size(mid,2)) * dy
+
+          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_yv_mask(mid,av), &
+                                  netcdf_data, start = (/ 1 /),    &
+                                  count = (/ mask_size(mid,2) /))
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          DEALLOCATE( netcdf_data )
+
+!
+!--       Write zu and zw data (vertical axes)
+          ALLOCATE( netcdf_data(mask_size(mid,3)) )
+
+          netcdf_data = zu( mask_k_global(mid,:mask_size(mid,3)) )
+
+          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
+                                  netcdf_data, start = (/ 1 /), &
+                                  count = (/ mask_size(mid,3) /) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          netcdf_data = zw( mask_k_global(mid,:mask_size(mid,3)) )
+
+          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
+                                  netcdf_data, start = (/ 1 /), &
+                                  count = (/ mask_size(mid,3) /) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          DEALLOCATE( netcdf_data )
+
+!
+!--       In case of non-flat topography write height information
+          IF ( TRIM( topography ) /= 'flat' )  THEN
+
+             ALLOCATE( netcdf_data_2d(mask_size(mid,1),mask_size(mid,2)) )
+             netcdf_data_2d = zu_s_inner( mask_i_global(mid,:mask_size(mid,1)),&
+                                          mask_j_global(mid,:mask_size(mid,2)) )
+
+             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),         &
+                                     id_var_zusi_mask(mid,av),    &
+                                     netcdf_data_2d,              &
+                                     start = (/ 1, 1 /),          &
+                                     count = (/ mask_size(mid,1), &
+                                                mask_size(mid,2) /) )
+             CALL handle_netcdf_error( 'netcdf', 9998 )
+
+             netcdf_data_2d = zw_w_inner( mask_i_global(mid,:mask_size(mid,1)),&
+                                          mask_j_global(mid,:mask_size(mid,2)) )
+
+             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),         &
+                                     id_var_zwwi_mask(mid,av),    &
+                                     netcdf_data_2d,              &
+                                     start = (/ 1, 1 /),          &
+                                     count = (/ mask_size(mid,1), &
+                                                mask_size(mid,2) /) )
+             CALL handle_netcdf_error( 'netcdf', 9998 )
+
+             DEALLOCATE( netcdf_data_2d )
+
+          ENDIF
+!
+!--       restore original parameter file_id (=formal parameter av) into av
+          av = file_id
+
+
+       CASE ( 'ma_ext' )
+
+!
+!--       decompose actual parameter file_id (=formal parameter av) into 
+!--       mid and av
+          file_id = av
+          IF ( file_id <= 140 )  THEN
+             mid = file_id - 120
+             av = 0
+          ELSE
+             mid = file_id - 140
+             av = 1
+          ENDIF
+
+!
+!--       Get the list of variables and compare with the actual run.
+!--       First var_list_old has to be reset, since GET_ATT does not assign
+!--       trailing blanks.
+          var_list_old = ' '
+          nc_stat = NF90_GET_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'VAR_LIST',&
+                                  var_list_old )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          var_list = ';'
+          i = 1
+          DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
+             var_list = TRIM(var_list) // TRIM( domask(mid,av,i) ) // ';'
+             i = i + 1
+          ENDDO
+
+          IF ( av == 0 )  THEN
+             var = '(mask)'
+          ELSE
+             var = '(mask_av)'
+          ENDIF
+
+          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
+             WRITE ( message_string, * ) 'NetCDF file for ', TRIM( var ), &
+                  ' data for mask', mid, ' from previous run found,', &
+                  '&but this file cannot be extended due to variable ', &
+                  'mismatch.&New file is created instead.'
+             CALL message( 'define_netcdf_header', 'PA9998', 0, 1, 0, 6, 0 )
+             extend = .FALSE.
+             RETURN
+          ENDIF
+
+!
+!--       Get and compare the number of vertical gridpoints
+          nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'zu', &
+                                    id_var_zu_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av),     &
+                                           id_var_zu_mask(mid,av),  &
+                                           dimids = id_dim_zu_mask_old )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+          id_dim_zu_mask(mid,av) = id_dim_zu_mask_old(1)
+
+          nc_stat = NF90_INQUIRE_DIMENSION( id_set_mask(mid,av),    &
+                                            id_dim_zu_mask(mid,av), &
+                                            len = nz_old )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          IF ( mask_size(mid,3) /= nz_old )  THEN
+             WRITE ( message_string, * ) 'NetCDF file for ', TRIM( var ), &
+                  ' data for mask', mid, ' from previous run found,', &
+                  '&but this file cannot be extended due to mismatch in ', &
+                  ' number of&vertical grid points.', &
+                  '&New file is created instead.'
+             CALL message( 'define_netcdf_header', 'PA9998', 0, 1, 0, 6, 0 )
+             extend = .FALSE.
+             RETURN
+          ENDIF
+
+!
+!--       Get the id of the time coordinate (unlimited coordinate) and its
+!--       last index on the file. The next time level is plmask..count+1.
+!--       The current time must be larger than the last output time
+!--       on the file.
+          nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'time', &
+                                    id_var_time_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av), &
+                                           id_var_time_mask(mid,av), &
+                                           dimids = id_dim_time_old )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+          id_dim_time_mask(mid,av) = id_dim_time_old(1)
+
+          nc_stat = NF90_INQUIRE_DIMENSION( id_set_mask(mid,av), &
+                                            id_dim_time_mask(mid,av), &
+                                            len = domask_time_count(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          nc_stat = NF90_GET_VAR( id_set_mask(mid,av), &
+                                  id_var_time_mask(mid,av), &
+                                  last_time_coordinate,              &
+                                  start = (/ domask_time_count(mid,av) /), &
+                                  count = (/ 1 /) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+
+          IF ( last_time_coordinate(1) >= simulated_time )  THEN
+             WRITE ( message_string, * ) 'NetCDF file for ', TRIM( var ), &
+                  ' data for mask', mid, ' from previous run found,', &
+                  '&but this file cannot be extended because the current ', &
+                  'output time&is less or equal than the last output time ', &
+                  'on this file.&New file is created instead.'
+             CALL message( 'define_netcdf_header', 'PA9998', 0, 1, 0, 6, 0 )
+             domask_time_count(mid,av) = 0
+             extend = .FALSE.
+             RETURN
+          ENDIF
+
+!
+!--       Dataset seems to be extendable.
+!--       Now get the variable ids.
+          i = 1
+          DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
+             nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), &
+                                       TRIM( domask(mid,av,i) ), &
+                                       id_var_domask(mid,av,i) )
+             CALL handle_netcdf_error( 'netcdf', 9998 )
+             i = i + 1
+          ENDDO
+
+!
+!--       Update the title attribute on file
+!--       In order to avoid 'data mode' errors if updated attributes are larger 
+!--       than their original size, NF90_PUT_ATT is called in 'define mode' 
+!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible 
+!--       performance loss due to data copying; an alternative strategy would be 
+!--       to ensure equal attribute size in a job chain. Maybe revise later.
+          IF ( av == 0 )  THEN
+             time_average_text = ' '
+          ELSE
+             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
+                                                            averaging_interval
+          ENDIF
+          nc_stat = NF90_REDEF( id_set_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'title', &
+                                  TRIM( run_description_header ) //    &
+                                  TRIM( time_average_text ) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+          nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
+          CALL handle_netcdf_error( 'netcdf', 9998 )
+          WRITE ( message_string, * ) 'NetCDF file for ', TRIM( var ), &
+               ' data for mask', mid, ' from previous run found.', &
+               '&This file will be extended.'
+          CALL message( 'define_netcdf_header', 'PA9998', 0, 0, 0, 6, 0 )
+!
+!--       restore original parameter file_id (=formal parameter av) into av
+          av = file_id
+
 
        CASE ( '3d_new' )
