Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 407)
+++ palm/trunk/SOURCE/Makefile	(revision 410)
@@ -5,4 +5,8 @@
 # -----------------
 # +set_particle_attributes, set_slicer_attributes_dvrp
+#
+# Branch revisions:
+# -----------------
+# +data_output_mask, init_masks, user_data_output_mask
 #
 # Former revisions:
@@ -55,14 +59,15 @@
 	close_file.f90 compute_vpt.f90 coriolis.f90 cpu_log.f90 \
 	cpu_statistics.f90 data_log.f90 data_output_dvrp.f90 \
-	data_output_profiles.f90 data_output_ptseries.f90 \
+	data_output_mask.f90 data_output_profiles.f90 data_output_ptseries.f90 \
 	data_output_spectra.f90 data_output_tseries.f90 data_output_2d.f90 \
 	data_output_3d.f90 diffusion_e.f90 diffusion_s.f90 diffusion_u.f90 \
 	diffusion_v.f90 diffusion_w.f90 diffusivities.f90 disturb_field.f90 \
-	disturb_heatflux.f90 eqn_state_seawater.f90 exchange_horiz.f90 exchange_horiz_2d.f90 \
+	disturb_heatflux.f90 eqn_state_seawater.f90 exchange_horiz.f90 \
+	exchange_horiz_2d.f90 \
 	fft_xy.f90 flow_statistics.f90 global_min_max.f90 \
 	header.f90 impact_of_latent_heat.f90 inflow_turbulence.f90 \
 	init_1d_model.f90 init_3d_model.f90 init_advec.f90 \
 	init_cloud_physics.f90 init_coupling.f90 init_dvrp.f90 init_grid.f90 \
-	init_ocean.f90 init_particles.f90 init_pegrid.f90 \
+	init_masks.f90 init_ocean.f90 init_particles.f90 init_pegrid.f90 \
 	init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \
 	interaction_droplets_ptq.f90 local_flush.f90 local_getenv.f90 \
@@ -83,5 +88,5 @@
 	user_check_data_output.f90 user_check_data_output_pr.f90 \
 	user_check_parameters.f90 user_data_output_2d.f90 \
-	user_data_output_3d.f90 user_data_output_dvrp.f90 \
+	user_data_output_3d.f90 user_data_output_dvrp.f90 user_data_output_mask.f90 \
 	user_define_netcdf_grid.f90 user_dvrp_coltab.f90 user_header.f90 \
 	user_init.f90 user_init_3d_model.f90 user_init_grid.f90 \
@@ -100,5 +105,6 @@
 	check_for_restart.o check_open.o check_parameters.o close_file.o \
 	compute_vpt.o coriolis.o cpu_log.o cpu_statistics.o data_log.o \
-	data_output_dvrp.o data_output_profiles.o data_output_ptseries.o \
+	data_output_dvrp.o data_output_mask.o data_output_profiles.o \
+	data_output_ptseries.o \
 	data_output_spectra.o data_output_tseries.o data_output_2d.o \
 	data_output_3d.o diffusion_e.o diffusion_s.o diffusion_u.o \
@@ -107,5 +113,5 @@
 	flow_statistics.o global_min_max.o header.o impact_of_latent_heat.o \
 	inflow_turbulence.o init_1d_model.o init_3d_model.o init_advec.o init_cloud_physics.o \
-	init_coupling.o init_dvrp.o init_grid.o init_ocean.o init_particles.o init_pegrid.o \
+	init_coupling.o init_dvrp.o init_grid.o init_masks.o init_ocean.o init_particles.o init_pegrid.o \
 	init_pt_anomaly.o init_rankine.o init_slope.o \
 	interaction_droplets_ptq.o local_flush.o local_getenv.o local_stop.o \
@@ -124,5 +130,5 @@
 	user_advec_particles.o user_check_data_output.o \
 	user_check_data_output_pr.o user_check_parameters.o \
-	user_data_output_2d.o user_data_output_3d.o user_data_output_dvrp.o \
+	user_data_output_2d.o user_data_output_3d.o user_data_output_mask.o user_data_output_dvrp.o \
 	user_define_netcdf_grid.o user_dvrp_coltab.o user_header.o \
 	user_init.o user_init_3d_model.o user_init_grid.o \
@@ -190,4 +196,5 @@
 data_log.o: modules.o
 data_output_dvrp.o: modules.o
+data_output_mask.o: modules.o
 data_output_profiles.o: modules.o
 data_output_ptseries.o: modules.o
@@ -220,4 +227,5 @@
 init_dvrp.o: modules.o
 init_grid.o: modules.o
+init_masks.o: modules.o
 init_ocean.o: modules.o eqn_state_seawater.o
 init_particles.o: modules.o random_function.o
@@ -282,4 +290,5 @@
 user_data_output_2d.o: modules.o user_module.o
 user_data_output_3d.o: modules.o user_module.o
+user_data_output_mask.o: modules.o user_module.o
 user_data_output_dvrp.o: modules.o user_module.o
 user_define_netcdf_grid.o: modules.o user_module.o
Index: palm/trunk/SOURCE/check_open.f90
===================================================================
--- palm/trunk/SOURCE/check_open.f90	(revision 407)
+++ palm/trunk/SOURCE/check_open.f90	(revision 410)
@@ -5,4 +5,8 @@
 ! -----------------
 ! 
+!
+! Branch revisions:
+! -----------------
+! masked data output
 !
 ! Former revisions:
@@ -60,5 +64,5 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=2)   ::  suffix
+    CHARACTER (LEN=2)   ::  mask_char, suffix
     CHARACTER (LEN=20)  ::  xtext = 'time in s'
     CHARACTER (LEN=30)  ::  filename
@@ -123,11 +127,15 @@
 
        CASE ( 15, 16, 17, 18, 19, 40:49, 50:59, 81:84, 101:107, 109, 111:113, &
-              116 )
+              116, 121:160 )
+	      
+	  IF ( .NOT. format_parallel_io) THEN
           
-          IF ( myid /= 0 )  THEN
-             WRITE( message_string, * ) 'opening file-id ',file_id, &
-                                        ' not allowed for PE ',myid
-             CALL message( 'check_open', 'PA0167', 2, 2, -1, 6, 1 )
-          ENDIF
+             IF ( myid /= 0 )  THEN
+                WRITE( message_string, * ) 'opening file-id ',file_id, &
+                                           ' not allowed for PE ',myid
+                CALL message( 'check_open', 'PA0167', 2, 2, -1, 6, 1 )
+             ENDIF
+	     
+	  ENDIF
 
        CASE ( 21, 22, 23 )
@@ -1233,7 +1241,117 @@
 
           ENDIF
+
+       CASE ( 121:160 )
+!
+!--       Set filename depending on unit number
+          IF ( file_id <= 140 )  THEN
+             mid = file_id - 120
+             WRITE ( mask_char,'(I2.2)')  mid
+             filename = 'DATA_MASK_' // mask_char // '_NETCDF' // coupling_char
+             av = 0
+          ELSE
+             mid = file_id - 140
+             WRITE ( mask_char,'(I2.2)')  mid
+             filename = 'DATA_MASK_' // mask_char // '_AV_NETCDF' // &
+                  coupling_char
+             av = 1
+          ENDIF
+!
+!--       Inquire, if there is a NetCDF file from a previuos run. This should
+!--       be opened for extension, if its dimensions and variables match the
+!--       actual run.
+          INQUIRE( FILE=filename, EXIST=netcdf_extend )
+
+          IF ( netcdf_extend )  THEN
+!
+!--          Open an existing NetCDF file for output 
+#if defined( __netcdf4 )
+             nc_stat = NF90_OPEN_PAR( filename, NF90_WRITE, comm2d, &
+                                      MPI_INFO_NULL, id_set_mask(mid,av) )
 #else
-
-       CASE ( 101:109, 111:113, 116 )
+             nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set_mask(mid,av) )
+#endif
+             CALL handle_netcdf_error( 'check_open', 9998 )
+!
+!--          Read header information and set all ids. If there is a mismatch
+!--          between the previuos and the actual run, netcdf_extend is returned
+!--          as .FALSE.
+             CALL define_netcdf_header( 'ma', netcdf_extend, file_id )
+
+!
+!--          Remove the local file, if it can not be extended
+             IF ( .NOT. netcdf_extend )  THEN
+                nc_stat = NF90_CLOSE( id_set_mask(mid,av) )
+                CALL handle_netcdf_error( 'check_open', 9998 )
+                CALL local_system('rm ' // TRIM( filename ) )
+             ENDIF
+
+          ENDIF          
+
+          IF ( .NOT. netcdf_extend )  THEN
+!
+!--          Create a new NetCDF output file
+             SELECT CASE ( nc_format_mask(mid,av) )
+
+                CASE ( 1 )
+                   nc_stat = NF90_CREATE( filename, NF90_NOCLOBBER, &
+                        id_set_mask(mid,av) )
+
+                CASE ( 2 )
+#if defined( __netcdf_64bit )
+                   nc_stat = NF90_CREATE( filename, OR( NF90_NOCLOBBER, &
+                        NF90_64BIT_OFFSET ), id_set_mask(mid,av) )
+#else
+                   WRITE( message_string, * ) 'NetCDF: no 64-bit ', &
+                        'offset format allowed on this machine'
+                   CALL message( 'check_open', 'PA9998', 0, 1, 0, 6, 0 )
+                   nc_stat = NF90_CREATE( filename, NF90_NOCLOBBER, &
+                        id_set_mask(mid,av) )
+#endif
+
+                CASE ( 3 )
+#if defined( __netcdf4 )
+                   nc_stat = NF90_CREATE_PAR( filename, OR( NF90_NOCLOBBER, &
+                        NF90_NETCDF4 ), comm2d, MPI_INFO_NULL, &
+                        id_set_mask(mid,av) )
+#else
+                   WRITE( message_string, * ) 'NetCDF: no NetCDF 4 ', &
+                        'format allowed on this machine'
+                   CALL message( 'check_open', 'PA9998', 0, 1, 0, 6, 0 )
+                   nc_stat = NF90_CREATE( filename, NF90_NOCLOBBER, &
+                        id_set_mask(mid,av) )
+#endif
+
+                CASE ( 4 )
+#if defined( __netcdf4 )
+                   nc_stat = NF90_CREATE_PAR( filename, OR( NF90_NOCLOBBER, &
+                      NF90_CLASSIC_MODEL ), comm2d, &
+                      MPI_INFO_NULL, id_set_mask(mid,av) )
+#else
+                   WRITE( message_string, * ) 'NetCDF: no NetCDF 4 (Classic ', &
+                        'model) format allowed on this machine'
+                   CALL message( 'check_open', 'PA9998', 0, 1, 0, 6, 0 )
+                   nc_stat = NF90_CREATE( filename, NF90_NOCLOBBER, &
+                        id_set_mask(mid,av) )
+#endif
+
+                CASE DEFAULT
+                   WRITE( message_string, * ) 'illegal NetCDF file format: ', &
+                        'nc_format_mask(mid=',mid,',av=',av,')=', &
+                        nc_format_mask(mid,av)
+                   CALL message( 'check_open', 'PA9998', 2, 2, 0, 6, 0 )
+
+             END SELECT
+             CALL handle_netcdf_error( 'check_open', 9998 )
+!
+!--          Define the header
+             CALL define_netcdf_header( 'ma', netcdf_extend, file_id )
+
+          ENDIF
+
+
+#else
+
+       CASE ( 101:109, 111:113, 116, 121:160 )
 
 !
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 407)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 410)
@@ -5,4 +5,8 @@
 ! -----------------
 ! Enabled passive scalar/humidity wall fluxes for non-flat topography
+!
+! Branch revisions:
+! -----------------
+! Masked data output
 !
 ! Former revisions:
@@ -1477,4 +1481,7 @@
        IF ( dt_do3d           == 9999999.9 )  dt_do3d           = dt_data_output
        IF ( dt_data_output_av == 9999999.9 )  dt_data_output_av = dt_data_output
+       DO  mid = 1, 20
+          IF ( dt_domask(mid) == 9999999.9 )  dt_domask(mid)    = dt_data_output
+       ENDDO
     ENDIF
 
@@ -1495,4 +1502,8 @@
     IF ( skip_time_data_output_av == 9999999.9 ) &
                                 skip_time_data_output_av = skip_time_data_output
+    DO  mid = 1, 20
+       IF ( skip_time_domask(mid) == 9999999.9 ) &
+                                skip_time_domask(mid)    = skip_time_data_output
+    ENDDO
 
 !
@@ -2647,4 +2658,26 @@
 
 !
+!-- Check mask conditions
+    IF ( masks < 0 .OR. masks > max_masks )  THEN
+       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ', &
+            '<= ', max_masks, ' (=max_masks)'
+       CALL message( 'check_parameters', 'PA9998', 1, 2, 0, 6, 0 )
+    ENDIF
+    IF ( masks > 0 )  THEN
+       mask_scale(1) = mask_scale_x
+       mask_scale(2) = mask_scale_y
+       mask_scale(3) = mask_scale_z
+       IF ( ANY( mask_scale <= 0.0 ) )  THEN
+          WRITE( message_string, * )  &
+               'illegal value: mask_scale_x, mask_scale_y and mask_scale_z', &
+               'must be > 0.0'
+          CALL message( 'check_parameters', 'PA9998', 1, 2, 0, 6, 0 )
+       ENDIF
+!
+!--    Generate masks for masked data output
+       CALL init_masks
+    ENDIF
+
+!
 !-- Check netcdf precison
     ldum = .FALSE.
Index: palm/trunk/SOURCE/close_file.f90
===================================================================
--- palm/trunk/SOURCE/close_file.f90	(revision 407)
+++ palm/trunk/SOURCE/close_file.f90	(revision 410)
@@ -5,4 +5,8 @@
 ! -----------------
 ! 
+!
+! Branch revisions:
+! -----------------
+! masked data output
 !
 ! Former revisions:
@@ -59,5 +63,5 @@
     CHARACTER (LEN=80)  ::  rtext, title, utext = '', xtext = '', ytext = ''
 
-    INTEGER ::  anzzeile, cranz, cross_count, cross_numbers, dimx, dimy, &
+    INTEGER ::  av, anzzeile, cranz, cross_count, cross_numbers, dimx, dimy, &
                 fid, file_id, j, k, legpos = 1, planz, timodex = 1
     INTEGER, DIMENSION(100) ::  klist, lstyle, cucol
@@ -94,5 +98,5 @@
 !
 !-- Close all open unit numbers
-    DO  fid = 1, 116
+    DO  fid = 1, 160
 
        IF ( openfile(fid)%opened .OR. openfile(fid)%opened_before )  THEN
@@ -493,4 +497,34 @@
                    nc_stat = NF90_CLOSE( id_set_3d(1) )
                    CALL handle_netcdf_error( 'close_file', 353 )
+                ENDIF
+
+             CASE ( 121:160 )
+	     
+	        IF ( format_parallel_io ) THEN
+!
+!--                decompose fid into mid and av
+                   IF ( fid <= 140 )  THEN
+                      mid = fid - 120
+                      av = 0
+                   ELSE
+                      mid = fid - 140
+                      av = 1
+                   ENDIF
+                   nc_stat = NF90_CLOSE( id_set_mask(mid,av) )
+                   CALL handle_netcdf_error( 'close_file', 9998 )
+		
+		ELSEIF ( myid == 0  .AND.  netcdf_output )  THEN
+!
+!--                decompose fid into mid and av
+                   IF ( fid <= 140 )  THEN
+                      mid = fid - 120
+                      av = 0
+                   ELSE
+                      mid = fid - 140
+                      av = 1
+                   ENDIF
+                   nc_stat = NF90_CLOSE( id_set_mask(mid,av) )
+                   CALL handle_netcdf_error( 'close_file', 9998 )
+		   
                 ENDIF
 
Index: palm/trunk/SOURCE/data_output_mask.f90
===================================================================
--- palm/trunk/SOURCE/data_output_mask.f90	(revision 410)
+++ palm/trunk/SOURCE/data_output_mask.f90	(revision 410)
@@ -0,0 +1,452 @@
+ SUBROUTINE data_output_mask( av )
+
+!------------------------------------------------------------------------------!
+! Actual revisions:
+! -----------------
+! Bugfix: to_be_resorted => s_av for time-averaged scalars
+! 
+!
+! Former revisions:
+! -----------------
+! $Id$
+!
+!
+! Description:
+! ------------
+! Masked data output in NetCDF format for current mask (current value of mid).
+!------------------------------------------------------------------------------!
+
+#if defined( __netcdf )
+    USE arrays_3d
+    USE averaging
+    USE cloud_parameters
+    USE control_parameters
+    USE cpulog
+    USE grid_variables
+    USE indices
+    USE interfaces
+    USE netcdf
+    USE netcdf_control
+    USE particle_attributes
+    USE pegrid
+
+    IMPLICIT NONE
+
+    INTEGER ::  av, ngp, file_id, i, if, is, j, k, l, n, psi, s, sender, &
+                ind(6)
+    LOGICAL ::  found, resorted
+    REAL    ::  mean_r, s_r3, s_r4
+    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  local_pf
+#if defined( __parallel )
+    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  total_pf
+#endif
+    REAL, DIMENSION(:,:,:), POINTER ::  to_be_resorted
+
+!
+!-- Return, if nothing to output
+    IF ( domask_no(mid,av) == 0 )  RETURN
+
+    CALL cpu_log (log_point(49),'data_output_mask','start')
+
+!!
+!-- Open output file.
+    IF ( format_parallel_io .AND. netcdf_output )  THEN
+       CALL check_open( 120+mid+av*max_masks )
+    ELSEIF (myid == 0  .AND.  netcdf_output ) THEN
+       CALL check_open( 120+mid+av*max_masks )
+    ENDIF 
+
+!
+!-- Allocate total and local output arrays.
+#if defined( __parallel )
+    IF ( myid == 0 )  THEN
+       ALLOCATE( total_pf(mask_size(mid,1),mask_size(mid,2),mask_size(mid,3)) )
+    ENDIF
+#endif
+    ALLOCATE( local_pf(mask_size_l(mid,1),mask_size_l(mid,2), &
+                       mask_size_l(mid,3)) )
+
+!
+!-- Update the NetCDF time axis.
+    domask_time_count(mid,av) = domask_time_count(mid,av) + 1
+    IF ( format_parallel_io .AND. netcdf_output )  THEN
+       nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_time_mask(mid,av), &
+                               (/ simulated_time /),                          &
+                               start = (/ domask_time_count(mid,av) /),       &
+                               count = (/ 1 /) )
+       CALL handle_netcdf_error( 'data_output_mask', 9998 )
+    ELSEIF ( myid == 0  .AND.  netcdf_output )  THEN
+       nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_time_mask(mid,av), &
+                               (/ simulated_time /),                          &
+                               start = (/ domask_time_count(mid,av) /),       &
+                               count = (/ 1 /) )
+       CALL handle_netcdf_error( 'data_output_mask', 9998 )
+    ENDIF
+
+!
+!-- Loop over all variables to be written.
+    if = 1
+
+    DO  WHILE ( domask(mid,av,if)(1:1) /= ' ' )
+!
+!--    Reallocate local_pf on PE 0 since its shape changes during MPI exchange
+       IF ( .NOT. format_parallel_io .AND. myid == 0 .AND. if > 1 )  THEN
+          DEALLOCATE( local_pf )
+          ALLOCATE( local_pf(mask_size_l(mid,1),mask_size_l(mid,2), &
+                             mask_size_l(mid,3)) )
+       ENDIF
+!
+!--    Store the variable chosen.
+       resorted = .FALSE.
+       SELECT CASE ( TRIM( domask(mid,av,if) ) )
+
+          CASE ( 'e' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => e
+             ELSE
+                to_be_resorted => e_av
+             ENDIF
+
+          CASE ( 'p' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => p
+             ELSE
+                to_be_resorted => p_av
+             ENDIF
+
+          CASE ( 'pc' )  ! particle concentration (requires ghostpoint exchange)
+             IF ( av == 0 )  THEN
+                tend = prt_count
+                CALL exchange_horiz( tend )
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) =  tend(mask_k(mid,k), &
+                                   mask_j(mid,j),mask_i(mid,i))
+                      ENDDO
+                   ENDDO
+                ENDDO
+                resorted = .TRUE.
+             ELSE
+                CALL exchange_horiz( pc_av )
+                to_be_resorted => pc_av
+             ENDIF
+
+          CASE ( 'pr' )  ! mean particle radius
+             IF ( av == 0 )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         psi = prt_start_index(k,j,i)
+                         s_r3 = 0.0
+                         s_r4 = 0.0
+                         DO  n = psi, psi+prt_count(k,j,i)-1
+                            s_r3 = s_r3 + particles(n)%radius**3
+                            s_r4 = s_r4 + particles(n)%radius**4
+                         ENDDO
+                         IF ( s_r3 /= 0.0 )  THEN
+                            mean_r = s_r4 / s_r3
+                         ELSE
+                            mean_r = 0.0
+                         ENDIF
+                         tend(k,j,i) = mean_r
+                      ENDDO
+                   ENDDO
+                ENDDO
+                CALL exchange_horiz( tend )
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) =  tend(mask_k(mid,k), &
+                                   mask_j(mid,j),mask_i(mid,i))
+                      ENDDO
+                   ENDDO
+                ENDDO
+                resorted = .TRUE.
+             ELSE
+                CALL exchange_horiz( pr_av )
+                to_be_resorted => pr_av
+             ENDIF
+
+          CASE ( 'pt' )
+             IF ( av == 0 )  THEN
+                IF ( .NOT. cloud_physics ) THEN
+                   to_be_resorted => pt
+                ELSE
+                   DO  i = 1, mask_size_l(mid,1)
+                      DO  j = 1, mask_size_l(mid,2)
+                         DO  k = 1, mask_size_l(mid,3)
+                            local_pf(i,j,k) =  &
+                                 pt(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) &
+                                 + l_d_cp * pt_d_t(mask_k(mid,k)) * &
+                                   ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                   resorted = .TRUE.
+                ENDIF
+             ELSE
+                to_be_resorted => pt_av
+             ENDIF
+
+          CASE ( 'q' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => q
+             ELSE
+                to_be_resorted => q_av
+             ENDIF
+             
+          CASE ( 'ql' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => ql
+             ELSE
+                to_be_resorted => ql_av
+             ENDIF
+
+          CASE ( 'ql_c' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => ql_c
+             ELSE
+                to_be_resorted => ql_c_av
+             ENDIF
+
+          CASE ( 'ql_v' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => ql_v
+             ELSE
+                to_be_resorted => ql_v_av
+             ENDIF
+
+          CASE ( 'ql_vp' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => ql_vp
+             ELSE
+                to_be_resorted => ql_vp_av
+             ENDIF
+
+          CASE ( 'qv' )
+             IF ( av == 0 )  THEN
+                DO  i = 1, mask_size_l(mid,1)
+                   DO  j = 1, mask_size_l(mid,2)
+                      DO  k = 1, mask_size_l(mid,3)
+                         local_pf(i,j,k) =  &
+                              q(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) -  &
+                              ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
+                      ENDDO
+                   ENDDO
+                ENDDO
+                resorted = .TRUE.
+             ELSE
+                to_be_resorted => qv_av
+             ENDIF
+
+          CASE ( 'rho' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rho
+             ELSE
+                to_be_resorted => rho_av
+             ENDIF
+             
+          CASE ( 's' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => q
+             ELSE
+                to_be_resorted => s_av
+             ENDIF
+             
+          CASE ( 'sa' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => sa
+             ELSE
+                to_be_resorted => sa_av
+             ENDIF
+             
+          CASE ( 'u' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => u
+             ELSE
+                to_be_resorted => u_av
+             ENDIF
+
+          CASE ( 'v' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => v
+             ELSE
+                to_be_resorted => v_av
+             ENDIF
+
+          CASE ( 'vpt' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => vpt
+             ELSE
+                to_be_resorted => vpt_av
+             ENDIF
+
+          CASE ( 'w' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => w
+             ELSE
+                to_be_resorted => w_av
+             ENDIF
+
+          CASE DEFAULT
+!
+!--          User defined quantity
+             CALL user_data_output_mask(av, domask(mid,av,if), found, local_pf )
+             resorted = .TRUE.
+
+             IF ( .NOT. found )  THEN
+                WRITE ( message_string, * ) 'no output available for: ', &
+                                            TRIM( domask(mid,av,if) )
+                CALL message( 'data_output_mask', 'PA9998', 0, 0, 0, 6, 0 )
+             ENDIF
+
+       END SELECT
+
+!
+!--    Resort the array to be output, if not done above
+       IF ( .NOT. resorted )  THEN
+          DO  i = 1, mask_size_l(mid,1)
+             DO  j = 1, mask_size_l(mid,2)
+                DO  k = 1, mask_size_l(mid,3)
+                   local_pf(i,j,k) =  to_be_resorted(mask_k(mid,k), &
+                                      mask_j(mid,j),mask_i(mid,i))
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
+
+!
+!--    I/O block. I/O methods are implemented
+!--    (1) for parallel execution
+!--     a. with NetCDF 4 parallel I/O-enabled library
+!--     b. with NetCDF 3 library
+!--    (2) for serial execution.
+!--    The choice of method depends on the correct setting of preprocessor
+!--    directives __parallel and __netcdf4 as well as on the parameter 
+!--    nc_format_mask(mid,av).
+#if defined( __parallel )
+#if defined( __netcdf4 )
+       IF ( nc_format_mask(mid,av) == 3 .OR. nc_format_mask(mid,av) == 4 )  THEN
+!
+!--       (1) a. Parallel I/O using NetCDF 4 (not yet tested)
+          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),  &
+               id_var_domask(mid,av,if),  &
+               local_pf,  &
+               start = (/ mask_start_l(mid,1), mask_start_l(mid,2),  &
+                          mask_start_l(mid,3), domask_time_count(mid,av) /),  &
+               count = (/ mask_size_l(mid,1), mask_size_l(mid,2),  &
+                          mask_size_l(mid,3), 1 /) )
+          CALL handle_netcdf_error( 'data_output_mask', 9998 )
+       ELSE
+#endif
+!
+!--       (1) b. Conventional I/O only through PE0
+!--       PE0 receives partial arrays from all processors of the respective mask
+!--       and outputs them. Here a barrier has to be set, because otherwise 
+!--       "-MPI- FATAL: Remote protocol queue full" may occur.
+          CALL MPI_BARRIER( comm2d, ierr )
+
+          ngp = mask_size_l(mid,1) * mask_size_l(mid,2) * mask_size_l(mid,3)
+          IF ( myid == 0 )  THEN
+!
+!--          Local array can be relocated directly.
+             total_pf( &
+               mask_start_l(mid,1):mask_start_l(mid,1)+mask_size_l(mid,1)-1, &
+               mask_start_l(mid,2):mask_start_l(mid,2)+mask_size_l(mid,2)-1, &
+               mask_start_l(mid,3):mask_start_l(mid,3)+mask_size_l(mid,3)-1 ) &
+               = local_pf
+!
+!--          Receive data from all other PEs.
+             DO  n = 1, numprocs-1
+!
+!--             Receive index limits first, then array.
+!--             Index limits are received in arbitrary order from the PEs.
+                CALL MPI_RECV( ind(1), 6, MPI_INTEGER, MPI_ANY_SOURCE, 0,  &
+                     comm2d, status, ierr )
+!
+!--             Not all PEs have data for the mask
+                IF ( ind(1) /= -9999 )  THEN
+                   ngp = ( ind(2)-ind(1)+1 ) * (ind(4)-ind(3)+1 ) *  &
+                         ( ind(6)-ind(5)+1 )
+                   sender = status(MPI_SOURCE)
+                   DEALLOCATE( local_pf )
+                   ALLOCATE(local_pf(ind(1):ind(2),ind(3):ind(4),ind(5):ind(6)))
+                   CALL MPI_RECV( local_pf(ind(1),ind(3),ind(5)), ngp,  &
+                        MPI_REAL, sender, 1, comm2d, status, ierr )
+                   total_pf(ind(1):ind(2),ind(3):ind(4),ind(5):ind(6)) &
+                        = local_pf
+                ENDIF
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),  &
+                  id_var_domask(mid,av,if), total_pf, &
+                  start = (/ 1, 1, 1, domask_time_count(mid,av) /), &
+                  count = (/ mask_size(mid,1), mask_size(mid,2), &
+                             mask_size(mid,3), 1 /) )
+             CALL handle_netcdf_error( 'data_output_mask', 9998 )
+
+          ELSE
+!
+!--          If at least part of the mask resides on the PE, send the index
+!--          limits for the target array, otherwise send -9999 to PE0.
+             IF ( mask_size_l(mid,1) > 0 .AND.  mask_size_l(mid,2) > 0 .AND. &
+                  mask_size_l(mid,3) > 0  ) &
+                  THEN
+                ind(1) = mask_start_l(mid,1)
+                ind(2) = mask_start_l(mid,1) + mask_size_l(mid,1) - 1
+                ind(3) = mask_start_l(mid,2)
+                ind(4) = mask_start_l(mid,2) + mask_size_l(mid,2) - 1
+                ind(5) = mask_start_l(mid,3)
+                ind(6) = mask_start_l(mid,3) + mask_size_l(mid,3) - 1
+             ELSE
+                ind(1) = -9999; ind(2) = -9999
+                ind(3) = -9999; ind(4) = -9999
+                ind(5) = -9999; ind(6) = -9999
+             ENDIF
+             CALL MPI_SEND( ind(1), 6, MPI_INTEGER, 0, 0, comm2d, ierr )
+!
+!--          If applicable, send data to PE0.
+             IF ( ind(1) /= -9999 )  THEN
+                CALL MPI_SEND( local_pf(1,1,1), ngp, MPI_REAL, 0, 1, comm2d, &
+                     ierr )
+             ENDIF
+          ENDIF
+!
+!--       A barrier has to be set, because otherwise some PEs may proceed too
+!--       fast so that PE0 may receive wrong data on tag 0.
+          CALL MPI_BARRIER( comm2d, ierr )
+#if defined( __netcdf4 )
+       ENDIF
+#endif
+#else
+!
+!--    (2) For serial execution of PALM, the single processor (PE0) holds all
+!--    data and writes them directly to file.
+       nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),  &
+            id_var_domask(mid,av,if),       &
+            local_pf(mask_size_l(mid,1),mask_size_l(mid,2), &
+                     mask_size_l(mid,3)),&
+            start = (/ 1, 1, 1, domask_time_count(mid,av) /), &
+            count = (/ mask_size_l(mid,1), mask_size_l(mid,2), &
+                       mask_size_l(mid,3), 1 /) )
+       CALL handle_netcdf_error( 'data_output_mask', 9998 )
+#endif
+
+       if = if + 1
+    ENDDO
+
+!
+!-- Deallocate temporary arrays.
+    DEALLOCATE( local_pf )
+#if defined( __parallel )
+    IF ( myid == 0 )  THEN
+       DEALLOCATE( total_pf )
+    ENDIF
+#endif
+
+
+    CALL cpu_log (log_point(49),'data_output_mask','stop','nobarrier')
+#endif
+
+ END SUBROUTINE data_output_mask
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 407)
+++ palm/trunk/SOURCE/header.f90	(revision 410)
@@ -6,4 +6,9 @@
 ! 
 !
+! Branch revisions:
+! -----------------
+! Masked data output: + dt_domask, mask_01~20_x|y|z, mask_01~20_x|y|z_loop, 
+! mask_scale|_x|y|z, masks, netcdf_format_mask[_av], skip_time_domask
+! 
 ! Former revisions:
 ! -----------------
@@ -105,11 +110,13 @@
     CHARACTER (LEN=70) ::  char1, char2, dopr_chr, &
                            do2d_xy, do2d_xz, do2d_yz, do3d_chr, &
-                           run_classification
+                           domask_chr, run_classification
     CHARACTER (LEN=86) ::  coordinates, gradients, learde, slices,  &
                            temperatures, ugcomponent, vgcomponent
     CHARACTER (LEN=85) ::  roben, runten
 
-    INTEGER ::  av, bh, blx, bly, bxl, bxr, byn, bys, ch, cwx, cwy, cxl, cxr, &
-                cyn, cys, i, ihost, io, j, l, ll, m, mpi_type
+    CHARACTER (LEN=1), DIMENSION(1:3) ::  dir = (/ 'x', 'y', 'z' /)
+
+    INTEGER ::  av, bh, blx, bly, bxl, bxr, byn, bys, ch, count, cwx, cwy,  &
+         cxl, cxr, cyn, cys, dim, i, ihost, io, j, l, ll, m, mpi_type
     REAL    ::  cpuseconds_per_simulated_second
 
@@ -988,4 +995,84 @@
        ENDIF
 
+    ENDDO
+
+!
+!-- masked arrays
+    IF ( masks > 0 )  WRITE ( io, 345 )  &
+         mask_scale_x, mask_scale_y, mask_scale_z
+    DO  mid = 1, masks
+       DO  av = 0, 1
+
+          i = 1
+          domask_chr = ''
+          DO  WHILE ( domask(mid,av,i) /= ' ' )
+             domask_chr = TRIM( domask_chr ) // ' ' //  &
+                          TRIM( domask(mid,av,i) ) // ','
+             i = i + 1
+          ENDDO
+
+          IF ( domask_chr /= '' )  THEN
+             IF ( av == 0 )  THEN
+                WRITE ( io, 346 )  '', mid
+             ELSE
+                WRITE ( io, 346 )  ' (time-averaged)', mid
+             ENDIF
+
+             output_format = ''
+             IF ( netcdf_output )  THEN
+                SELECT CASE ( nc_format_mask(mid,av) )
+                   CASE ( 1 )
+                      output_format = 'netcdf (classic format)'
+                   CASE ( 2 )
+                      output_format = 'netcdf (64bit offset format)'
+                   CASE ( 3 )
+                      output_format = 'netcdf (NetCDF 4 format)'
+                   CASE ( 4 )
+                      output_format = 'netcdf (NetCDF 4 classic model format)'
+                END SELECT
+             ENDIF
+             WRITE ( io, 344 )  output_format
+
+             IF ( av == 0 )  THEN
+                WRITE ( io, 347 )  domask_chr, dt_domask(mid)
+             ELSE
+                WRITE ( io, 348 )  domask_chr, dt_data_output_av, &
+                                   averaging_interval, dt_averaging_input
+             ENDIF
+
+             IF ( av == 0 )  THEN
+                IF ( skip_time_domask(mid) /= 0.0 )  THEN
+                   WRITE ( io, 339 )  skip_time_domask(mid)
+                ENDIF
+             ELSE
+                IF ( skip_time_data_output_av /= 0.0 )  THEN
+                   WRITE ( io, 339 )  skip_time_data_output_av
+                ENDIF
+             ENDIF
+!
+!--          output locations
+             DO  dim = 1, 3
+                IF ( mask(mid,dim,1) >= 0.0 )  THEN
+                   count = 0
+                   DO  WHILE ( mask(mid,dim,count+1) >= 0.0 )
+                      count = count + 1
+                   ENDDO
+                   WRITE ( io, 349 )  dir(dim), dir(dim), mid, dir(dim), &
+                                      mask(mid,dim,:count)
+                ELSEIF ( mask_loop(mid,dim,1) < 0.0 .AND.  &
+                         mask_loop(mid,dim,2) < 0.0 .AND.  &
+                         mask_loop(mid,dim,3) == 0.0 )  THEN
+                   WRITE ( io, 350 )  dir(dim), dir(dim)
+                ELSEIF ( mask_loop(mid,dim,3) == 0.0 )  THEN
+                   WRITE ( io, 351 )  dir(dim), dir(dim), mid, dir(dim), &
+                                      mask_loop(mid,dim,1:2)
+                ELSE
+                   WRITE ( io, 351 )  dir(dim), dir(dim), mid, dir(dim), &
+                                      mask_loop(mid,dim,1:3)
+                ENDIF
+             ENDDO
+          ENDIF
+
+       ENDDO
     ENDDO
 
@@ -1619,4 +1706,23 @@
             '       Upper output limit at    ',F8.2,' m  (GP ',I4,')'/)
 344 FORMAT ('       Output format: ',A/)
+345 FORMAT (/'    Scaling lengths for output locations of all subsequent mask IDs:',/ &
+            '       mask_scale_x (in x-direction): ',F9.3, ' m',/ &
+            '       mask_scale_y (in y-direction): ',F9.3, ' m',/ &
+            '       mask_scale_z (in z-direction): ',F9.3, ' m' )
+346 FORMAT (/'    Masked data output',A,' for mask ID ',I2, ':')
+347 FORMAT ('       Variables: ',A/ &
+            '       Output every             ',F8.2,' s')
+348 FORMAT ('       Variables: ',A/ &
+            '       Output every             ',F8.2,' s'/ &
+            '       Time averaged over       ',F8.2,' s'/ &
+            '       Averaging input every    ',F8.2,' s')
+349 FORMAT (/'       Output locations in ',A,'-direction in multiples of ', &
+            'mask_scale_',A,' predefined by array mask_',I2.2,'_',A,':'/ &
+            13('       ',8(F8.2,',')/) )
+350 FORMAT (/'       Output locations in ',A,'-direction: ', &
+            'all gridpoints along ',A,'-direction (default).' )
+351 FORMAT (/'       Output locations in ',A,'-direction in multiples of ', &
+            'mask_scale_',A,' constructed from array mask_',I2.2,'_',A,'_loop:'/ &
+            '          loop begin:',F8.2,', end:',F8.2,', stride:',F8.2 )
 #if defined( __dvrp_graphics )
 360 FORMAT ('    Plot-Sequence with dvrp-software:'/ &
Index: palm/trunk/SOURCE/init_masks.f90
===================================================================
--- palm/trunk/SOURCE/init_masks.f90	(revision 410)
+++ palm/trunk/SOURCE/init_masks.f90	(revision 410)
@@ -0,0 +1,650 @@
+ SUBROUTINE init_masks
+
+!------------------------------------------------------------------------------!
+! Actual revisions:
+! -----------------
+! 
+!
+! To do:
+! ------
+! 
+!
+! Former revisions:
+! -----------------
+! $Id$
+!
+!
+! Description:
+! ------------
+! Initialize masked data output
+!------------------------------------------------------------------------------!
+
+    USE arrays_3d
+    USE control_parameters
+    USE grid_variables
+    USE indices
+    USE netcdf_control
+    USE particle_attributes
+    USE pegrid
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=6)   ::  var
+    CHARACTER (LEN=7)   ::  unit
+    CHARACTER (LEN=10), DIMENSION(max_masks,100) ::  do_mask, do_mask_user
+
+    INTEGER :: i, ilen, ind(6), ind_array(1), j, k, n, sender
+    INTEGER, DIMENSION(:), ALLOCATABLE ::  tmp_array
+
+    LOGICAL ::  found
+!
+!-- Allocation and initialization
+    ALLOCATE( tmp_array( MAX(nx,ny,nz)+1 ) )
+
+    ALLOCATE( mask_i(max_masks,nxr-nxl+1), &
+              mask_j(max_masks,nyn-nys+1), &
+              mask_k(max_masks,nzt-nzb+1) )
+!
+!-- internal mask arrays ("mask,dimension,selection")
+    ALLOCATE( mask(max_masks,3,mask_xyz_dimension), &
+              mask_loop(max_masks,3,3) )
+!
+!-- Store netcdf file formats for masked data output in one shared array
+    nc_format_mask(:,0) = netcdf_format_mask
+    nc_format_mask(:,1) = netcdf_format_mask_av
+    
+    IF ( ANY( nc_format_mask == 3 ) .OR. ANY( nc_format_mask == 4 )) THEN
+       format_parallel_io = .TRUE.
+    ENDIF
+
+    IF ( ANY( nc_format_mask < 1 ) .OR. ANY( nc_format_mask > 4 ) )  THEN
+       WRITE( message_string, * )  'illegal value: netcdf_format_mask and ', &
+            'netcdf_format_mask_av must be either 1, 2, 3 or 4'
+       IF ( ANY( nc_format_mask == 3 ) .OR. ANY( nc_format_mask == 4 ) )  THEN
+          message_string = TRIM( message_string ) // '&NetCDF file formats '// &
+               '3 (NetCDF 4) and 4 (NetCDF 4 Classic model)'// &
+               '&are currently not supported (not yet tested).'
+       ENDIF
+       CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+    ENDIF
+!
+!-- Store data output parameters for masked data output in few shared arrays
+    do_mask( 1,:) = data_output_mask_01; do_mask( 2,:) = data_output_mask_02
+    do_mask( 3,:) = data_output_mask_03; do_mask( 4,:) = data_output_mask_04
+    do_mask( 5,:) = data_output_mask_05; do_mask( 6,:) = data_output_mask_06
+    do_mask( 7,:) = data_output_mask_07; do_mask( 8,:) = data_output_mask_08
+    do_mask( 9,:) = data_output_mask_09; do_mask(10,:) = data_output_mask_10
+    do_mask(11,:) = data_output_mask_11; do_mask(12,:) = data_output_mask_12
+    do_mask(13,:) = data_output_mask_13; do_mask(14,:) = data_output_mask_14
+    do_mask(15,:) = data_output_mask_15; do_mask(16,:) = data_output_mask_16
+    do_mask(17,:) = data_output_mask_17; do_mask(18,:) = data_output_mask_18
+    do_mask(19,:) = data_output_mask_19; do_mask(20,:) = data_output_mask_20
+    do_mask_user( 1,:) = data_output_mask_user_01
+    do_mask_user( 2,:) = data_output_mask_user_02
+    do_mask_user( 3,:) = data_output_mask_user_03
+    do_mask_user( 4,:) = data_output_mask_user_04
+    do_mask_user( 5,:) = data_output_mask_user_05
+    do_mask_user( 6,:) = data_output_mask_user_06
+    do_mask_user( 7,:) = data_output_mask_user_07
+    do_mask_user( 8,:) = data_output_mask_user_08
+    do_mask_user( 9,:) = data_output_mask_user_09
+    do_mask_user(10,:) = data_output_mask_user_10
+    do_mask_user(11,:) = data_output_mask_user_11
+    do_mask_user(12,:) = data_output_mask_user_12
+    do_mask_user(13,:) = data_output_mask_user_13
+    do_mask_user(14,:) = data_output_mask_user_14
+    do_mask_user(15,:) = data_output_mask_user_15
+    do_mask_user(16,:) = data_output_mask_user_16
+    do_mask_user(17,:) = data_output_mask_user_17
+    do_mask_user(18,:) = data_output_mask_user_18
+    do_mask_user(19,:) = data_output_mask_user_19
+    do_mask_user(20,:) = data_output_mask_user_20
+    mask( 1,1,:) = mask_01_x; mask( 2,1,:) = mask_02_x
+    mask( 3,1,:) = mask_03_x; mask( 4,1,:) = mask_04_x
+    mask( 5,1,:) = mask_05_x; mask( 6,1,:) = mask_06_x
+    mask( 7,1,:) = mask_07_x; mask( 8,1,:) = mask_08_x
+    mask( 9,1,:) = mask_09_x; mask(10,1,:) = mask_10_x
+    mask(11,1,:) = mask_11_x; mask(12,1,:) = mask_12_x
+    mask(13,1,:) = mask_13_x; mask(14,1,:) = mask_14_x
+    mask(15,1,:) = mask_15_x; mask(16,1,:) = mask_16_x
+    mask(17,1,:) = mask_17_x; mask(18,1,:) = mask_18_x
+    mask(19,1,:) = mask_19_x; mask(20,1,:) = mask_20_x
+    mask( 1,2,:) = mask_01_y; mask( 2,2,:) = mask_02_y
+    mask( 3,2,:) = mask_03_y; mask( 4,2,:) = mask_04_y
+    mask( 5,2,:) = mask_05_y; mask( 6,2,:) = mask_06_y
+    mask( 7,2,:) = mask_07_y; mask( 8,2,:) = mask_08_y
+    mask( 9,2,:) = mask_09_y; mask(10,2,:) = mask_10_y
+    mask(11,2,:) = mask_11_y; mask(12,2,:) = mask_12_y
+    mask(13,2,:) = mask_13_y; mask(14,2,:) = mask_14_y
+    mask(15,2,:) = mask_15_y; mask(16,2,:) = mask_16_y
+    mask(17,2,:) = mask_17_y; mask(18,2,:) = mask_18_y
+    mask(19,2,:) = mask_19_y; mask(20,2,:) = mask_20_y
+    mask( 1,3,:) = mask_01_z; mask( 2,3,:) = mask_02_z
+    mask( 3,3,:) = mask_03_z; mask( 4,3,:) = mask_04_z
+    mask( 5,3,:) = mask_05_z; mask( 6,3,:) = mask_06_z
+    mask( 7,3,:) = mask_07_z; mask( 8,3,:) = mask_08_z
+    mask( 9,3,:) = mask_09_z; mask(10,3,:) = mask_10_z
+    mask(11,3,:) = mask_11_z; mask(12,3,:) = mask_12_z
+    mask(13,3,:) = mask_13_z; mask(14,3,:) = mask_14_z
+    mask(15,3,:) = mask_15_z; mask(16,3,:) = mask_16_z
+    mask(17,3,:) = mask_17_z; mask(18,3,:) = mask_18_z
+    mask(19,3,:) = mask_19_z; mask(20,3,:) = mask_20_z
+    mask_loop( 1,1,:) = mask_01_x_loop; mask_loop( 2,1,:) = mask_02_x_loop
+    mask_loop( 3,1,:) = mask_03_x_loop; mask_loop( 4,1,:) = mask_04_x_loop
+    mask_loop( 5,1,:) = mask_05_x_loop; mask_loop( 6,1,:) = mask_06_x_loop
+    mask_loop( 7,1,:) = mask_07_x_loop; mask_loop( 8,1,:) = mask_08_x_loop
+    mask_loop( 9,1,:) = mask_09_x_loop; mask_loop(10,1,:) = mask_10_x_loop
+    mask_loop(11,1,:) = mask_11_x_loop; mask_loop(12,1,:) = mask_12_x_loop
+    mask_loop(13,1,:) = mask_13_x_loop; mask_loop(14,1,:) = mask_14_x_loop
+    mask_loop(15,1,:) = mask_15_x_loop; mask_loop(16,1,:) = mask_16_x_loop
+    mask_loop(17,1,:) = mask_17_x_loop; mask_loop(18,1,:) = mask_18_x_loop
+    mask_loop(19,1,:) = mask_19_x_loop; mask_loop(20,1,:) = mask_20_x_loop
+    mask_loop( 1,2,:) = mask_01_y_loop; mask_loop( 2,2,:) = mask_02_y_loop
+    mask_loop( 3,2,:) = mask_03_y_loop; mask_loop( 4,2,:) = mask_04_y_loop
+    mask_loop( 5,2,:) = mask_05_y_loop; mask_loop( 6,2,:) = mask_06_y_loop
+    mask_loop( 7,2,:) = mask_07_y_loop; mask_loop( 8,2,:) = mask_08_y_loop
+    mask_loop( 9,2,:) = mask_09_y_loop; mask_loop(10,2,:) = mask_10_y_loop
+    mask_loop(11,2,:) = mask_11_y_loop; mask_loop(12,2,:) = mask_12_y_loop
+    mask_loop(13,2,:) = mask_13_y_loop; mask_loop(14,2,:) = mask_14_y_loop
+    mask_loop(15,2,:) = mask_15_y_loop; mask_loop(16,2,:) = mask_16_y_loop
+    mask_loop(17,2,:) = mask_17_y_loop; mask_loop(18,2,:) = mask_18_y_loop
+    mask_loop(19,2,:) = mask_19_y_loop; mask_loop(20,2,:) = mask_20_y_loop
+    mask_loop( 1,3,:) = mask_01_z_loop; mask_loop( 2,3,:) = mask_02_z_loop
+    mask_loop( 3,3,:) = mask_03_z_loop; mask_loop( 4,3,:) = mask_04_z_loop
+    mask_loop( 5,3,:) = mask_05_z_loop; mask_loop( 6,3,:) = mask_06_z_loop
+    mask_loop( 7,3,:) = mask_07_z_loop; mask_loop( 8,3,:) = mask_08_z_loop
+    mask_loop( 9,3,:) = mask_09_z_loop; mask_loop(10,3,:) = mask_10_z_loop
+    mask_loop(11,3,:) = mask_11_z_loop; mask_loop(12,3,:) = mask_12_z_loop
+    mask_loop(13,3,:) = mask_13_z_loop; mask_loop(14,3,:) = mask_14_z_loop
+    mask_loop(15,3,:) = mask_15_z_loop; mask_loop(16,3,:) = mask_16_z_loop
+    mask_loop(17,3,:) = mask_17_z_loop; mask_loop(18,3,:) = mask_18_z_loop
+    mask_loop(19,3,:) = mask_19_z_loop; mask_loop(20,3,:) = mask_20_z_loop
+    mask_i = -1; mask_j = -1; mask_k = -1
+!
+!-- Global arrays are required by define_netcdf_header.
+    IF ( format_parallel_io )  THEN
+       ALLOCATE( mask_i_global(max_masks,nx+1), &
+                 mask_j_global(max_masks,ny+1), &
+                 mask_k_global(max_masks,nz+1) )
+       mask_i_global = -1; mask_j_global = -1; mask_k_global = -1
+    ELSEIF ( myid == 0 )  THEN
+       ALLOCATE( mask_i_global(max_masks,nx+1), &
+                 mask_j_global(max_masks,ny+1), &
+                 mask_k_global(max_masks,nz+1) )
+       mask_i_global = -1; mask_j_global = -1; mask_k_global = -1
+    ENDIF
+
+!
+!-- Determine variable names for each mask
+    DO  mid = 1, masks
+!
+!--    Append user-defined data output variables to the standard data output
+       IF ( do_mask_user(mid,1) /= ' ' )  THEN
+          i = 1
+          DO  WHILE ( do_mask(mid,i) /= ' '  .AND.  i <= 100 )
+             i = i + 1
+          ENDDO
+          j = 1
+          DO  WHILE ( do_mask_user(mid,j) /= ' '  .AND.  j <= 100 )
+             IF ( i > 100 )  THEN
+                WRITE ( message_string, * ) 'number of output quantitities ', &
+                     'given by data_output_mask and data_output_mask_user ', &
+                     'exceeds the limit of 100'
+                CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+             ENDIF
+             do_mask(mid,i) = do_mask_user(mid,j)
+             i = i + 1
+             j = j + 1
+          ENDDO
+       ENDIF
+
+!
+!--    Check and set steering parameters for mask data output and averaging
+       i   = 1
+       DO  WHILE ( do_mask(mid,i) /= ' '  .AND.  i <= 100 )
+!
+!--       Check for data averaging
+          ilen = LEN_TRIM( do_mask(mid,i) )
+          j = 0                                              ! no data averaging
+          IF ( ilen > 3 )  THEN
+             IF ( do_mask(mid,i)(ilen-2:ilen) == '_av' )  THEN
+                j = 1                                           ! data averaging
+                do_mask(mid,i) = do_mask(mid,i)(1:ilen-3)
+             ENDIF
+          ENDIF
+          var = TRIM( do_mask(mid,i) )
+!
+!--       Check for allowed value and set units
+          SELECT CASE ( TRIM( var ) )
+
+             CASE ( 'e' )
+                IF ( constant_diffusion )  THEN
+                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
+                        '" requires constant_diffusion = .FALSE.'
+                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'm2/s2'
+
+             CASE ( 'pc', 'pr' )
+                IF ( .NOT. particle_advection )  THEN
+                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
+                        '" requires a "particles_par"-NAMELIST in the ', &
+                        'parameter file (PARIN)'
+                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+                ENDIF
+                IF ( TRIM( var ) == 'pc' )  unit = 'number'
+                IF ( TRIM( var ) == 'pr' )  unit = 'm'
+
+             CASE ( 'q', 'vpt' )
+                IF ( .NOT. humidity )  THEN
+                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
+                        '" requires humidity = .TRUE.'
+                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+                ENDIF
+                IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
+                IF ( TRIM( var ) == 'vpt' )  unit = 'K'
+
+             CASE ( 'ql' )
+                IF ( .NOT. ( cloud_physics  .OR.  cloud_droplets ) )  THEN
+                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
+                        '" requires cloud_physics = .TRUE. or cloud_droplets', &
+                        ' = .TRUE.'
+                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'kg/kg'
+
+             CASE ( 'ql_c', 'ql_v', 'ql_vp' )
+                IF ( .NOT. cloud_droplets )  THEN
+                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
+                        '" requires cloud_droplets = .TRUE.'
+                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+                ENDIF
+                IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
+                IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
+                IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
+
+             CASE ( 'qv' )
+                IF ( .NOT. cloud_physics )  THEN
+                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
+                        '" requires cloud_physics = .TRUE.'
+                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'kg/kg'
+
+             CASE ( 'rho' )
+                IF ( .NOT. ocean )  THEN
+                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
+                        '" requires ocean = .TRUE.'
+                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'kg/m3'
+
+             CASE ( 's' )
+                IF ( .NOT. passive_scalar )  THEN
+                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
+                        '" requires passive_scalar = .TRUE.'
+                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'conc'
+
+             CASE ( 'sa' )
+                IF ( .NOT. ocean )  THEN
+                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
+                        '" requires ocean = .TRUE.'
+                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'psu'
+
+             CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*' )
+                WRITE ( message_string, * ) 'illegal value for data_',&
+                     'output: "', TRIM( var ), '" is only allowed', &
+                     'for horizontal cross section'
+                CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+!                IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
+!                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
+!                         '" requires cloud_physics = .TRUE.'
+!                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+!                ENDIF
+!                IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
+!                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
+!                         '" requires precipitation = .TRUE.'
+!                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+!                ENDIF
+!                IF ( TRIM( var ) == 'pra*'  .AND.  j == 1 )  THEN
+!                   WRITE ( message_string, * ) 'temporal averaging of ', &
+!                         ' precipitation amount "', TRIM( var ),         &
+!                         '" not possible' 
+!                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+!                ENDIF
+!                IF ( TRIM( var ) == 'prr*'  .AND.  .NOT. precipitation )  THEN
+!                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
+!                         '" requires precipitation = .TRUE.'
+!                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+!                ENDIF
+!
+!                IF ( TRIM( var ) == 'u*'   )  unit = 'm/s'
+!                IF ( TRIM( var ) == 't*'   )  unit = 'K'
+!                IF ( TRIM( var ) == 'lwp*' )  unit = 'kg/kg*m'
+!                IF ( TRIM( var ) == 'pra*' )  unit = 'mm'
+!                IF ( TRIM( var ) == 'prr*' )  unit = 'mm/s'
+!                IF ( TRIM( var ) == 'z0*'  )  unit = 'm'
+
+             CASE ( 'p', 'pt', 'u', 'v', 'w' )
+                IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
+                IF ( TRIM( var ) == 'pt' )  unit = 'K'
+                IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
+                IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
+                IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
+                CONTINUE
+
+             CASE DEFAULT
+                CALL user_check_data_output( var, unit )
+
+                IF ( unit == 'illegal' )  THEN
+                   IF ( do_mask_user(mid,1) /= ' ' )  THEN
+                      WRITE ( message_string, * ) 'illegal value for data_',&
+                           'output or data_output_user: "', &
+                           TRIM( do_mask(mid,i) ), '"'
+                      CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+                   ELSE
+                      WRITE ( message_string, * ) 'illegal value for data_',&
+                           'output: "', TRIM( do_mask(mid,i) ), '"'
+                      CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+                   ENDIF
+                ENDIF
+
+          END SELECT
+!
+!--       Set the internal steering parameters appropriately
+          domask_no(mid,j)                    = domask_no(mid,j) + 1
+          domask(mid,j,domask_no(mid,j))      = do_mask(mid,i)
+          domask_unit(mid,j,domask_no(mid,j)) = unit
+
+          IF ( j == 1 )  THEN
+!
+!--          Check, if variable is already subject to averaging
+             found = .FALSE.
+             DO  k = 1, doav_n
+                IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
+             ENDDO
+
+             IF ( .NOT. found )  THEN
+                doav_n = doav_n + 1
+                doav(doav_n) = var
+             ENDIF
+          ENDIF
+
+          i = i + 1
+
+       ENDDO   ! do_mask(mid,i)
+    ENDDO   ! mid
+
+
+!
+!-- Determine mask locations for each mask
+    DO  mid = 1, masks
+!
+!--    Set local masks for each subdomain along all three dimensions
+       CALL set_mask_locations( 1, dx, 'dx', nx, 'nx', nxl, nxr )
+       CALL set_mask_locations( 2, dy, 'dy', ny, 'ny', nys, nyn )
+       CALL set_mask_locations( 3, dz, 'dz', nz, 'nz', nzb, nzt )
+!
+!--    Set global masks along all three dimensions (required by
+!--    define_netcdf_header).
+#if defined( __parallel )
+!
+!--    PE0 receives partial arrays from all processors of the respective mask
+!--    and outputs them. Here a barrier has to be set, because otherwise 
+!--    "-MPI- FATAL: Remote protocol queue full" may occur.
+       CALL MPI_BARRIER( comm2d, ierr )
+
+       IF ( myid == 0 )  THEN
+!
+!--       Local arrays can be relocated directly.
+          mask_i_global(mid,mask_start_l(mid,1): &
+                       mask_start_l(mid,1)+mask_size_l(mid,1)-1) = &
+                       mask_i(mid,:mask_size_l(mid,1))
+          mask_j_global(mid,mask_start_l(mid,2): &
+                       mask_start_l(mid,2)+mask_size_l(mid,2)-1) = &
+                       mask_j(mid,:mask_size_l(mid,2))
+          mask_k_global(mid,mask_start_l(mid,3): &
+                       mask_start_l(mid,3)+mask_size_l(mid,3)-1) = &
+                       mask_k(mid,:mask_size_l(mid,3))
+!
+!--       Receive data from all other PEs.
+          DO  n = 1, numprocs-1
+!
+!--          Receive index limits first, then arrays.
+!--          Index limits are received in arbitrary order from the PEs.
+             CALL MPI_RECV( ind(1), 6, MPI_INTEGER, MPI_ANY_SOURCE, 0,  &
+                  comm2d, status, ierr )
+!
+!--          Not all PEs have data for the mask.
+             IF ( ind(1) /= -9999 )  THEN
+                sender = status(MPI_SOURCE)
+                CALL MPI_RECV( tmp_array(ind(1)), ind(2)-ind(1)+1,  &
+                               MPI_INTEGER, sender, 1, comm2d, status, ierr )
+                mask_i_global(mid,ind(1):ind(2)) = tmp_array(ind(1):ind(2))
+                CALL MPI_RECV( tmp_array(ind(3)), ind(4)-ind(3)+1,  &
+                               MPI_INTEGER, sender, 2, comm2d, status, ierr )
+                mask_j_global(mid,ind(3):ind(4)) = tmp_array(ind(3):ind(4))
+                CALL MPI_RECV( tmp_array(ind(5)), ind(6)-ind(5)+1,  &
+                               MPI_INTEGER, sender, 3, comm2d, status, ierr )
+                mask_k_global(mid,ind(5):ind(6)) = tmp_array(ind(5):ind(6))
+             ENDIF
+          ENDDO
+
+       ELSE
+!
+!--       If at least part of the mask resides on the PE, send the index limits
+!--       for the target array, otherwise send -9999 to PE0.
+          IF ( mask_size_l(mid,1) > 0 .AND.  mask_size_l(mid,2) > 0 .AND. &
+               mask_size_l(mid,3) > 0  ) &
+               THEN
+             ind(1) = mask_start_l(mid,1)
+             ind(2) = mask_start_l(mid,1) + mask_size_l(mid,1) - 1
+             ind(3) = mask_start_l(mid,2)
+             ind(4) = mask_start_l(mid,2) + mask_size_l(mid,2) - 1
+             ind(5) = mask_start_l(mid,3)
+             ind(6) = mask_start_l(mid,3) + mask_size_l(mid,3) - 1
+          ELSE
+             ind(1) = -9999; ind(2) = -9999
+             ind(3) = -9999; ind(4) = -9999
+             ind(5) = -9999; ind(6) = -9999
+          ENDIF
+          CALL MPI_SEND( ind(1), 6, MPI_INTEGER, 0, 0, comm2d, ierr )
+!
+!--       If applicable, send data to PE0.
+          IF ( ind(1) /= -9999 )  THEN
+             tmp_array(:mask_size_l(mid,1)) = mask_i(mid,:mask_size_l(mid,1))
+             CALL MPI_SEND( tmp_array(1), mask_size_l(mid,1),  &
+                            MPI_INTEGER, 0, 1, comm2d, ierr )
+             tmp_array(:mask_size_l(mid,2)) = mask_j(mid,:mask_size_l(mid,2))
+             CALL MPI_SEND( tmp_array(1), mask_size_l(mid,2),  &
+                            MPI_INTEGER, 0, 2, comm2d, ierr )
+             tmp_array(:mask_size_l(mid,3)) = mask_k(mid,:mask_size_l(mid,3))
+             CALL MPI_SEND( tmp_array(1), mask_size_l(mid,3),  &
+                            MPI_INTEGER, 0, 3, comm2d, ierr )
+          ENDIF
+       ENDIF
+!
+!--    A barrier has to be set, because otherwise some PEs may proceed too fast
+!--    so that PE0 may receive wrong data on tag 0.
+       CALL MPI_BARRIER( comm2d, ierr )
+       
+       IF ( format_parallel_io )  THEN  
+         
+          CALL MPI_BCAST( mask_i_global(mid,:), nx+1, MPI_INTEGER, 0, comm2d, ierr )
+          CALL MPI_BCAST( mask_j_global(mid,:), ny+1, MPI_INTEGER, 0, comm2d, ierr )
+          CALL MPI_BCAST( mask_k_global(mid,:), nz+1, MPI_INTEGER, 0, comm2d, ierr ) 
+	  
+       ENDIF
+
+#else
+!
+!--    Local arrays can be relocated directly.
+       mask_i_global(mid,:) = mask_i(mid,:)
+       mask_j_global(mid,:) = mask_j(mid,:)
+       mask_k_global(mid,:) = mask_k(mid,:)
+#endif
+    ENDDO   ! mid
+
+    DEALLOCATE( tmp_array )
+!
+!-- Internal mask arrays cannot be deallocated on PE 0 because they are
+!-- required for header output on PE 0.
+    IF ( myid /= 0 )  DEALLOCATE( mask, mask_loop )
+
+ CONTAINS
+    SUBROUTINE set_mask_locations( dim, dxyz, dxyz_string, nxyz, nxyz_string, &
+                                   lb, ub )
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+! Set local mask for each subdomain along 'dim' direction.
+!------------------------------------------------------------------------------!
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN=2) :: dxyz_string, nxyz_string
+       INTEGER  ::  count, count_l, dim, m, loop_begin, loop_end, loop_stride, &
+                    lb, nxyz, ub
+       REAL     ::  dxyz, ddxyz, tmp1, tmp2
+
+       count = 0;  count_l = 0;  ddxyz = 1.0 / dxyz;  tmp1 = 0.0;  tmp2 = 0.0
+
+       IF ( mask(mid,dim,1) >= 0.0 )  THEN
+!
+!--       use predefined mask_* array
+          DO  WHILE ( mask(mid,dim,count+1) >= 0.0 )
+             count = count + 1
+             IF ( dim == 1 .OR. dim == 2 )  THEN
+                m = NINT( mask(mid,dim,count) * mask_scale(dim) * ddxyz )
+             ELSEIF ( dim == 3 )  THEN
+                ind_array =  &
+                     MINLOC( ABS( mask(mid,dim,count) * mask_scale(dim) - zw ) )
+                m = ind_array(1) - 1 + nzb  ! MINLOC uses lower array bound 1
+             ENDIF
+             IF ( m > nxyz )  THEN
+                WRITE ( message_string, '(I3,A,I3,A,I1,3A,I3)' )  &
+                     m,' in mask ',mid,' along dimension ',dim,  &
+                     ' exceeds ',nxyz_string,' = ',nxyz
+                CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( m >= lb .AND. m <= ub )  THEN
+                IF ( count_l == 0 )  mask_start_l(mid,dim) = count
+                count_l = count_l + 1
+                IF ( dim == 1 )  THEN
+                   mask_i(mid,count_l) = m
+                ELSEIF ( dim == 2 )  THEN
+                   mask_j(mid,count_l) = m
+                ELSEIF ( dim == 3 )  THEN
+                   mask_k(mid,count_l) = m
+                ENDIF
+             ENDIF
+             IF ( count == mask_xyz_dimension )  EXIT
+          ENDDO
+          mask_size(mid,dim)   = count
+          mask_size_l(mid,dim) = count_l
+
+       ELSE
+!
+!--       use predefined mask_loop_* array, or use the default (all grid points
+!--       along this direction)
+          IF ( mask_loop(mid,dim,1) < 0.0 )  THEN
+             tmp1 = mask_loop(mid,dim,1)
+             mask_loop(mid,dim,1) = 0.0   ! (default)
+          ENDIF
+          IF ( dim == 1 .OR. dim == 2 )  THEN
+             IF ( mask_loop(mid,dim,2) < 0.0 )  THEN
+                tmp2 = mask_loop(mid,dim,2)
+                mask_loop(mid,dim,2) = nxyz*dxyz / mask_scale(dim)   ! (default)
+             ENDIF
+             IF ( MAXVAL( mask_loop(mid,dim,1:2) )  &
+                  > nxyz * dxyz / mask_scale(dim) )  THEN
+                WRITE ( message_string, '(2(A,I3,A,I1,A,F9.3),5A,I1,A,F9.3)' ) &
+                     'mask_loop(',mid,',',dim,',1)=',mask_loop(mid,dim,1), &
+                     ' and/or mask_loop(',mid,',',dim,',2)=', &
+                     mask_loop(mid,dim,2),'&exceed ', &
+                     nxyz_string,'*',dxyz_string,'/mask_scale(',dim,')=', &
+                     nxyz*dxyz/mask_scale(dim)
+                CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+             ENDIF
+             loop_begin  = NINT( mask_loop(mid,dim,1) * mask_scale(dim) * ddxyz )
+             loop_end    = NINT( mask_loop(mid,dim,2) * mask_scale(dim) * ddxyz )
+             loop_stride = NINT( mask_loop(mid,dim,3) * mask_scale(dim) * ddxyz )
+          ELSEIF ( dim == 3 )  THEN
+             IF ( mask_loop(mid,dim,2) < 0.0 )  THEN
+                tmp2 = mask_loop(mid,dim,2)
+                mask_loop(mid,dim,2) = zw(nz) / mask_scale(dim)   ! (default)
+             ENDIF
+             IF ( MAXVAL( mask_loop(mid,dim,1:2) )  &
+                  > zw(nz) / mask_scale(dim) )  THEN
+                WRITE ( message_string, '(2(A,I3,A,I1,A,F9.3),A,I1,A,F9.3)' ) &
+                     'mask_loop(',mid,',',dim,',1)=',mask_loop(mid,dim,1), &
+                     ' and/or mask_loop(',mid,',',dim,',2)=', &
+                     mask_loop(mid,dim,2),'&exceed zw(nz)/mask_scale(',dim, &
+                     ')=',zw(nz)/mask_scale(dim)
+                CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
+             ENDIF
+             ind_array =  &
+                  MINLOC( ABS( mask_loop(mid,dim,1) * mask_scale(dim) - zw ) )
+             loop_begin =  &
+                  ind_array(1) - 1 + nzb ! MINLOC uses lower array bound 1
+             ind_array =  &
+                  MINLOC( ABS( mask_loop(mid,dim,2) * mask_scale(dim) - zw ) )
+             loop_end = ind_array(1) - 1 + nzb ! MINLOC uses lower array bound 1
+!
+!--          The following line assumes a constant vertical grid spacing within
+!--          the vertical mask range; it fails for vertical grid stretching.
+!--          Maybe revise later. Issue warning but continue execution.
+             loop_stride = NINT( mask_loop(mid,dim,3) * mask_scale(dim) * ddxyz )
+
+             IF ( mask_loop(mid,dim,2) * mask_scale(dim) > dz_stretch_level )  &
+                  THEN
+                WRITE ( message_string, '(A,I3,A,I1,A,F9.3,A,F8.2,3A)' )  &
+                     'mask_loop(',mid,dim,',2)=', mask_loop(mid,dim,2),&
+                     ' exceeds dz_stretch_level=',dz_stretch_level, &
+                     '.&Vertical mask locations will not ', &
+                     'match the desired heights&within the stretching ', &
+                     'region. Recommendation: use mask instead of mask_loop.'
+                CALL message( 'init_masks', 'PA9998', 0, 1, 0, 6, 0 )
+             ENDIF
+
+          ENDIF
+!
+!--       If necessary, reset mask_loop(mid,dim,1) and mask_loop(mid,dim,2).
+          IF ( tmp1 < 0.0 )  mask_loop(mid,dim,1) = tmp1
+          IF ( tmp2 < 0.0 )  mask_loop(mid,dim,2) = tmp2
+!
+!--       The default stride +/-1 (every grid point) applies if 
+!--       mask_loop(mid,dim,3) is not specified (its default is zero).
+          IF ( loop_stride == 0 )  THEN
+             IF ( loop_end >= loop_begin )  THEN
+                loop_stride =  1
+             ELSE
+                loop_stride = -1
+             ENDIF
+          ENDIF
+          DO  m = loop_begin, loop_end, loop_stride
+             count = count + 1
+             IF ( m >= lb .AND. m <= ub )  THEN
+                IF ( count_l == 0 )  mask_start_l(mid,dim) = count
+                count_l = count_l + 1
+                IF ( dim == 1 )  THEN
+                   mask_i(mid,count_l) = m
+                ELSEIF ( dim == 2 )  THEN
+                   mask_j(mid,count_l) = m
+                ELSEIF ( dim == 3 )  THEN
+                   mask_k(mid,count_l) = m
+                ENDIF
+             ENDIF
+          ENDDO
+          mask_size(mid,dim)   = count
+          mask_size_l(mid,dim) = count_l
+       ENDIF
+
+    END SUBROUTINE set_mask_locations
+
+ END SUBROUTINE init_masks
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 407)
+++ palm/trunk/SOURCE/modules.f90	(revision 410)
@@ -8,4 +8,9 @@
 ! initial data assignments to some dvrp arrays changed due to error messages
 ! from gfortran compiler
+!
+! Branch revisions:
+! -----------------
+! Masked data output: + dt_domask, mask_01~20_x|y|z, mask_01~20_x|y|z_loop, 
+! mask_scale|_x|y|z, masks, netcdf_format_mask[_av], skip_time_domask
 !
 ! Former revisions:
@@ -281,4 +286,6 @@
     TYPE(file_status), DIMENSION(200) :: openfile = file_status(.FALSE.,.FALSE.)
 
+
+    INTEGER, PARAMETER :: mask_xyz_dimension = 100, max_masks = 20
 
     CHARACTER (LEN=1)    ::  cycle_mg = 'w', timestep_reason = ' '
@@ -322,8 +329,31 @@
     CHARACTER (LEN=11), DIMENSION(100) ::  data_output = ' ',    &
                                            data_output_user = ' ', doav = ' '
+    CHARACTER (LEN=10), DIMENSION(100) ::  &
+         data_output_mask_01      = ' ', data_output_mask_02      = ' ',  &
+         data_output_mask_03      = ' ', data_output_mask_04      = ' ',  &
+         data_output_mask_05      = ' ', data_output_mask_06      = ' ',  &
+         data_output_mask_07      = ' ', data_output_mask_08      = ' ',  &
+         data_output_mask_09      = ' ', data_output_mask_10      = ' ',  &
+         data_output_mask_11      = ' ', data_output_mask_12      = ' ',  &
+         data_output_mask_13      = ' ', data_output_mask_14      = ' ',  &
+         data_output_mask_15      = ' ', data_output_mask_16      = ' ',  &
+         data_output_mask_17      = ' ', data_output_mask_18      = ' ',  &
+         data_output_mask_19      = ' ', data_output_mask_20      = ' ',  &
+         data_output_mask_user_01 = ' ', data_output_mask_user_02 = ' ',  &
+         data_output_mask_user_03 = ' ', data_output_mask_user_04 = ' ',  &
+         data_output_mask_user_05 = ' ', data_output_mask_user_06 = ' ',  &
+         data_output_mask_user_07 = ' ', data_output_mask_user_08 = ' ',  &
+         data_output_mask_user_09 = ' ', data_output_mask_user_10 = ' ',  &
+         data_output_mask_user_11 = ' ', data_output_mask_user_12 = ' ',  &
+         data_output_mask_user_13 = ' ', data_output_mask_user_14 = ' ',  &
+         data_output_mask_user_15 = ' ', data_output_mask_user_16 = ' ',  &
+         data_output_mask_user_17 = ' ', data_output_mask_user_18 = ' ',  &
+         data_output_mask_user_19 = ' ', data_output_mask_user_20 = ' '
+
     CHARACTER (LEN=10), DIMENSION(300) ::  data_output_pr = ' '
     CHARACTER (LEN=10), DIMENSION(200) ::  data_output_pr_user = ' '
-    CHARACTER (LEN=20), DIMENSION(10)  ::  netcdf_precision = ' '
-
+    CHARACTER (LEN=20), DIMENSION(50)  ::  netcdf_precision = ' '
+
+    CHARACTER (LEN=10), DIMENSION(max_masks,0:1,100) ::  domask = ' '
     CHARACTER (LEN=10), DIMENSION(0:1,100) ::  do2d = ' ', do3d = ' '
 
@@ -341,7 +371,8 @@
                 inflow_disturbance_begin = -1, inflow_disturbance_end = -1, &
                 intermediate_timestep_count, intermediate_timestep_count_max, &
-                iran = -1234567, last_dt_change = 0, maximum_grid_level, &
+                iran = -1234567, last_dt_change = 0, masks = 0, &
+                maximum_grid_level, &
                 max_pr_user = 0, mgcycles = 0, mg_cycles = -1, &
-                mg_switch_to_pe0_level = 0, ngsrb = 2, nsor = 20, &
+                mg_switch_to_pe0_level = 0, mid, ngsrb = 2, nsor = 20, &
                 nsor_ini = 100, n_sor, normalizing_region = 0, &
                 nz_do1d, nz_do3d = -9999, outflow_damping_width = -1, &
@@ -354,5 +385,11 @@
                 do2d_xz_time_count(0:1), do2d_yz_time_count(0:1), &
                 do3d_no(0:1) = 0, do3d_time_count(0:1), &
+                domask_no(max_masks,0:1) = 0, domask_time_count(max_masks,0:1),&
                 lad_vertical_gradient_level_ind(10) = -9999, &
+                mask_size(max_masks,3) = -1, mask_size_l(max_masks,3) = -1, &
+                mask_start_l(max_masks,3) = -1, &
+                nc_format_mask(1:max_masks,0:1), &
+                netcdf_format_mask(1:max_masks) = 1, &
+                netcdf_format_mask_av(1:max_masks) = 1, &
                 pt_vertical_gradient_level_ind(10) = -9999, &
                 q_vertical_gradient_level_ind(10) = -9999, &
@@ -364,4 +401,8 @@
 
     INTEGER, DIMENSION(:), ALLOCATABLE ::  grid_level_count
+
+    INTEGER, DIMENSION(:,:), ALLOCATABLE   ::  mask_i, mask_j, mask_k
+    INTEGER, DIMENSION(:,:), ALLOCATABLE   ::  &
+                mask_i_global, mask_j_global, mask_k_global
 
     LOGICAL ::  adjust_mixing_length = .FALSE., avs_output = .FALSE., &
@@ -381,5 +422,6 @@
                 disturbance_created = .FALSE., &
                 first_call_advec_particles = .TRUE., &
-                force_print_header = .FALSE., galilei_transformation = .FALSE.,&
+                force_print_header = .FALSE., format_parallel_io = .FALSE., &
+		galilei_transformation = .FALSE.,&
                 humidity = .FALSE., humidity_remote = .FALSE., &
                 inflow_l = .FALSE., inflow_n = .FALSE., &
@@ -439,5 +481,7 @@
              km_damp_max = -1.0, lad_surface = 0.0,  &
              leaf_surface_concentration = 0.0, long_filter_factor = 0.0, &
-             maximum_cpu_time_allowed = 0.0, molecular_viscosity = 1.461E-5, &
+             mask_scale_x = 1.0, mask_scale_y = 1.0, mask_scale_z = 1.0, &
+             maximum_cpu_time_allowed = 0.0,  &
+             molecular_viscosity = 1.461E-5, &
              old_dt = 1.0E-10, omega = 7.29212E-5, omega_sor = 1.8, &
              overshoot_limit_e = 0.0, overshoot_limit_pt = 0.0, &
@@ -480,6 +524,7 @@
     REAL ::  do2d_xy_last_time(0:1) = -1.0, do2d_xz_last_time(0:1) = -1.0, &
              do2d_yz_last_time(0:1) = -1.0, dpdxy(1:2) = 0.0, &
-             lad_vertical_gradient(10) = 0.0, &
+             dt_domask(max_masks) = 9999999.9, lad_vertical_gradient(10) = 0.0,&
              lad_vertical_gradient_level(10) = -9999999.9, &
+             mask_scale(3), &
              pt_vertical_gradient(10) = 0.0, &
              pt_vertical_gradient_level(10) = -9999999.9, &
@@ -489,5 +534,7 @@
              s_vertical_gradient_level(10) = -1.0, &
              sa_vertical_gradient(10) = 0.0, &
-             sa_vertical_gradient_level(10) = -9999999.9, threshold(20) = 0.0, &
+             sa_vertical_gradient_level(10) = -9999999.9, &
+             skip_time_domask(max_masks) = 9999999.9, threshold(20) = 0.0, &
+             time_domask(max_masks) = 0.0, &
              tsc(10) = (/ 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /), &
              ug_vertical_gradient(10) = 0.0, &
@@ -502,4 +549,84 @@
     REAL, DIMENSION(:), ALLOCATABLE ::  dp_smooth_factor
 
+    REAL, DIMENSION(mask_xyz_dimension) ::  &
+       mask_01_x = -1.0, mask_02_x = -1.0, mask_03_x = -1.0, mask_04_x = -1.0, &
+       mask_05_x = -1.0, mask_06_x = -1.0, mask_07_x = -1.0, mask_08_x = -1.0, &
+       mask_09_x = -1.0, mask_10_x = -1.0, mask_11_x = -1.0, mask_12_x = -1.0, &
+       mask_13_x = -1.0, mask_14_x = -1.0, mask_15_x = -1.0, mask_16_x = -1.0, &
+       mask_17_x = -1.0, mask_18_x = -1.0, mask_19_x = -1.0, mask_20_x = -1.0, &
+       mask_01_y = -1.0, mask_02_y = -1.0, mask_03_y = -1.0, mask_04_y = -1.0, &
+       mask_05_y = -1.0, mask_06_y = -1.0, mask_07_y = -1.0, mask_08_y = -1.0, &
+       mask_09_y = -1.0, mask_10_y = -1.0, mask_11_y = -1.0, mask_12_y = -1.0, &
+       mask_13_y = -1.0, mask_14_y = -1.0, mask_15_y = -1.0, mask_16_y = -1.0, &
+       mask_17_y = -1.0, mask_18_y = -1.0, mask_19_y = -1.0, mask_20_y = -1.0, &
+       mask_01_z = -1.0, mask_02_z = -1.0, mask_03_z = -1.0, mask_04_z = -1.0, &
+       mask_05_z = -1.0, mask_06_z = -1.0, mask_07_z = -1.0, mask_08_z = -1.0, &
+       mask_09_z = -1.0, mask_10_z = -1.0, mask_11_z = -1.0, mask_12_z = -1.0, &
+       mask_13_z = -1.0, mask_14_z = -1.0, mask_15_z = -1.0, mask_16_z = -1.0, &
+       mask_17_z = -1.0, mask_18_z = -1.0, mask_19_z = -1.0, mask_20_z = -1.0
+    REAL, DIMENSION(3) ::  &
+       mask_01_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_02_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_03_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_04_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_05_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_06_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_07_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_08_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_09_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_10_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_11_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_12_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_13_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_14_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_15_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_16_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_17_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_18_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_19_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_20_x_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_01_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_02_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_03_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_04_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_05_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_06_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_07_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_08_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_09_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_10_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_11_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_12_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_13_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_14_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_15_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_16_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_17_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_18_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_19_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_20_y_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_01_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_02_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_03_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_04_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_05_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_06_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_07_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_08_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_09_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_10_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_11_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_12_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_13_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_14_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_15_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_16_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_17_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_18_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_19_z_loop = (/ -1.0, -1.0, 0.0 /), &
+       mask_20_z_loop = (/ -1.0, -1.0, 0.0 /)
+!
+!--    internal mask arrays ("mask,dimension,selection")
+       REAL, DIMENSION(:,:,:), ALLOCATABLE ::  mask, mask_loop
 
     SAVE
@@ -832,4 +959,5 @@
 !------------------------------------------------------------------------------!
 
+    USE control_parameters
 #if defined( __netcdf )
     USE netcdf
@@ -934,6 +1062,7 @@
                 id_var_zw_xy, id_var_zw_xz, id_var_zw_yz, id_var_zw_3d
 
-    INTEGER, DIMENSION(10)  ::  id_var_dospx, id_var_dospy, nc_precision 
+    INTEGER, DIMENSION(10)  ::  id_var_dospx, id_var_dospy
     INTEGER, DIMENSION(20)  ::  id_var_prt
+    INTEGER, DIMENSION(50)  ::  nc_precision
     INTEGER, DIMENSION(dopr_norm_num) ::  id_var_norm_dopr
 
@@ -944,5 +1073,19 @@
     INTEGER, DIMENSION(dots_max,0:9)   ::  id_var_dots
 
+!
+!-- masked output
+    CHARACTER (LEN=7), DIMENSION(max_masks,0:1,100) ::  domask_unit
+
     LOGICAL ::  output_for_t0 = .FALSE.
+
+    INTEGER, DIMENSION(1:max_masks,0:1) ::  id_dim_time_mask, id_dim_x_mask, &
+                id_dim_xu_mask, id_dim_y_mask, id_dim_yv_mask, id_dim_zu_mask, &
+                id_dim_zw_mask, &
+                id_set_mask, &
+                id_var_time_mask, id_var_x_mask, id_var_xu_mask, &
+                id_var_y_mask, id_var_yv_mask, id_var_zu_mask, id_var_zw_mask, &
+                id_var_zusi_mask, id_var_zwwi_mask
+
+    INTEGER, DIMENSION(1:max_masks,0:1,100)         ::  id_var_domask
 
     SAVE
Index: palm/trunk/SOURCE/netcdf.f90
===================================================================
--- palm/trunk/SOURCE/netcdf.f90	(revision 407)
+++ 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' )
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 407)
+++ palm/trunk/SOURCE/parin.f90	(revision 410)
@@ -5,4 +5,9 @@
 ! -----------------
 ! +wall_humidityflux, wall_scalarflux
+!
+! Branch revisions:
+! -----------------
+! Masked data output: + dt_domask, mask_01~20_x|y|z, mask_01~20_x|y|z_loop, 
+! mask_scale_x|y|z, masks, netcdf_format_mask[_av], skip_time_domask
 !
 ! Former revisions:
@@ -108,5 +113,5 @@
              lad_vertical_gradient, lad_vertical_gradient_level, &
              leaf_surface_concentration, long_filter_factor, &
-             loop_optimization, mixing_length_1d, momentum_advec, &
+             loop_optimization, mixing_length_1d, masks, momentum_advec, &
              netcdf_precision, nsor_ini, nx, ny, &
              nz, ocean, omega, outflow_damping_width, overshoot_limit_e, &
@@ -140,5 +145,16 @@
                        cross_normalized_y, cross_profiles, cross_ts_uymax, &
                        cross_ts_uymin, cross_xtext, cycle_mg, data_output, &
-                       data_output_format, data_output_pr, &
+                       data_output_format, &
+                       data_output_mask_01, data_output_mask_02, &
+                       data_output_mask_03, data_output_mask_04, &
+                       data_output_mask_05, data_output_mask_06, &
+                       data_output_mask_07, data_output_mask_08, &
+                       data_output_mask_09, data_output_mask_10, &
+                       data_output_mask_11, data_output_mask_12, &
+                       data_output_mask_13, data_output_mask_14, &
+                       data_output_mask_15, data_output_mask_16, &
+                       data_output_mask_17, data_output_mask_18, &
+                       data_output_mask_19, data_output_mask_20, &
+                       data_output_pr, &
                        data_output_2d_on_each_pe, disturbance_amplitude, &
                        disturbance_energy_limit, disturbance_level_b, &
@@ -146,9 +162,48 @@
                        do3d_compress, do3d_comp_prec, dt, dt_averaging_input, &
                        dt_averaging_input_pr, dt_coupling, dt_data_output, &
-                       dt_data_output_av, dt_disturb, dt_dopr, &
+                       dt_data_output_av, dt_disturb, dt_domask, dt_dopr, &
                        dt_dopr_listing, dt_dots, dt_do2d_xy, dt_do2d_xz, &
                        dt_do2d_yz, dt_do3d, dt_max, dt_restart, dt_run_control,&
-                       end_time, force_print_header, mg_cycles, &
-                       mg_switch_to_pe0_level, netcdf_64bit, netcdf_64bit_3d, &
+                       end_time, force_print_header, &
+                       mask_scale_x, mask_scale_y, mask_scale_z, &
+                       mask_01_x, mask_02_x, mask_03_x, mask_04_x, &
+                       mask_05_x, mask_06_x, mask_07_x, mask_08_x, &
+                       mask_09_x, mask_10_x, mask_11_x, mask_12_x, &
+                       mask_13_x, mask_14_x, mask_15_x, mask_16_x, &
+                       mask_17_x, mask_18_x, mask_19_x, mask_20_x, &
+                       mask_01_y, mask_02_y, mask_03_y, mask_04_y, &
+                       mask_05_y, mask_06_y, mask_07_y, mask_08_y, &
+                       mask_09_y, mask_10_y, mask_11_y, mask_12_y, &
+                       mask_13_y, mask_14_y, mask_15_y, mask_16_y, &
+                       mask_17_y, mask_18_y, mask_19_y, mask_20_y, &
+                       mask_01_z, mask_02_z, mask_03_z, mask_04_z, &
+                       mask_05_z, mask_06_z, mask_07_z, mask_08_z, &
+                       mask_09_z, mask_10_z, mask_11_z, mask_12_z, &
+                       mask_13_z, mask_14_z, mask_15_z, mask_16_z, &
+                       mask_17_z, mask_18_z, mask_19_z, mask_20_z, &
+                       mask_01_x_loop, mask_02_x_loop, mask_03_x_loop, &
+                       mask_04_x_loop, mask_05_x_loop, mask_06_x_loop, &
+                       mask_07_x_loop, mask_08_x_loop, mask_09_x_loop, &
+                       mask_10_x_loop, mask_11_x_loop, mask_12_x_loop, &
+                       mask_13_x_loop, mask_14_x_loop, mask_15_x_loop, &
+                       mask_16_x_loop, mask_17_x_loop, mask_18_x_loop, &
+                       mask_19_x_loop, mask_20_x_loop, &
+                       mask_01_y_loop, mask_02_y_loop, mask_03_y_loop, &
+                       mask_04_y_loop, mask_05_y_loop, mask_06_y_loop, &
+                       mask_07_y_loop, mask_08_y_loop, mask_09_y_loop, &
+                       mask_10_y_loop, mask_11_y_loop, mask_12_y_loop, &
+                       mask_13_y_loop, mask_14_y_loop, mask_15_y_loop, &
+                       mask_16_y_loop, mask_17_y_loop, mask_18_y_loop, &
+                       mask_19_y_loop, mask_20_y_loop, &
+                       mask_01_z_loop, mask_02_z_loop, mask_03_z_loop, &
+                       mask_04_z_loop, mask_05_z_loop, mask_06_z_loop, &
+                       mask_07_z_loop, mask_08_z_loop, mask_09_z_loop, &
+                       mask_10_z_loop, mask_11_z_loop, mask_12_z_loop, &
+                       mask_13_z_loop, mask_14_z_loop, mask_15_z_loop, &
+                       mask_16_z_loop, mask_17_z_loop, mask_18_z_loop, &
+                       mask_19_z_loop, mask_20_z_loop, &
+                       mg_cycles, mg_switch_to_pe0_level, &
+                       netcdf_64bit, netcdf_64bit_3d, &
+                       netcdf_format_mask, netcdf_format_mask_av, &
                        ngsrb, normalizing_region, npex, npey, nsor, nz_do3d, &
                        omega_sor, &
@@ -160,5 +215,5 @@
                        skip_time_data_output_av, skip_time_dopr, &
                        skip_time_do2d_xy, skip_time_do2d_xz, &
-                       skip_time_do2d_yz, skip_time_do3d, &
+                       skip_time_do2d_yz, skip_time_do3d, skip_time_domask, &
                        termination_time_needed, use_prior_plot1d_parameters, &
                        z_max_do1d, z_max_do1d_normalized, z_max_do2d
Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 407)
+++ palm/trunk/SOURCE/production_e.f90	(revision 410)
@@ -635,5 +635,5 @@
        USE indices
        USE statistics
-
+USE pegrid
        IMPLICIT NONE
 
@@ -988,4 +988,7 @@
 
                 DO  k = nzb_diff_s_inner(j,i), nzt_diff
+IF ( pt(k,j,i) == 0.0 .OR. dd2zu(k) == 0.0 )  &
+WRITE (*,'(6(I4,1X),F6.2,1X,F6.2)')  &
+myid, i, j, nzb_diff_s_inner(j,i), nzt_diff, k, pt(k,j,i), dd2zu(k)
                    tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
                                        ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
Index: palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/read_3d_binary.f90	(revision 407)
+++ palm/trunk/SOURCE/read_3d_binary.f90	(revision 410)
@@ -181,5 +181,5 @@
           DO  k = 1, j
              WRITE (9,*) 'k = ', k
-             WRITE (9,'(6(A,I2))')  'nxlfa = ', nxlfa(files_to_be_opened,k), &
+             WRITE (9,'(6(A,I4))')  'nxlfa = ', nxlfa(files_to_be_opened,k), &
                          ' nxrfa = ', nxrfa(files_to_be_opened,k), &
                          ' offset_xa = ', offset_xa(files_to_be_opened,k), &
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 407)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 410)
@@ -5,4 +5,8 @@
 ! -----------------
 ! Bugfix: exchange of ghost points for prho included
+!
+! Branch revisions:
+! -----------------
+! masked data output
 !
 ! Former revisions:
@@ -324,4 +328,9 @@
           time_do3d          = time_do3d        + dt_3d
        ENDIF
+       DO  mid = 1, masks
+          IF ( simulated_time >= skip_time_domask(mid) )  THEN
+             time_domask(mid)= time_domask(mid) + dt_3d
+          ENDIF
+       ENDDO
        time_dvrp          = time_dvrp        + dt_3d
        IF ( simulated_time >= skip_time_dosp )  THEN
@@ -399,5 +408,5 @@
 
 !
-!--    Sum-up 3d-arrays for later output of time-averaged data
+!--    Sum-up 3d-arrays for later output of time-averaged 2d/3d/masked data
        IF ( averaging_interval /= 0.0  .AND.                                &
             ( dt_data_output_av - time_do_av ) <= averaging_interval  .AND. &
@@ -495,5 +504,15 @@
 
 !
-!--    Output of time-averaged 2d/3d-data
+!--    masked data output
+       DO  mid = 1, masks
+          IF ( time_domask(mid) >= dt_domask(mid) )  THEN
+             CALL data_output_mask( 0 )
+             time_domask(mid) = MOD( time_domask(mid),  &
+                                     MAX( dt_domask(mid), dt_3d ) )
+          ENDIF
+       ENDDO
+
+!
+!--    Output of time-averaged 2d/3d/masked data
        IF ( time_do_av >= dt_data_output_av )  THEN
           CALL average_3d_data
@@ -502,4 +521,7 @@
           CALL data_output_2d( 'yz', 1 )
           CALL data_output_3d( 1 )
+          DO  mid = 1, masks
+             CALL data_output_mask( 1 )
+          ENDDO
           time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
        ENDIF
Index: palm/trunk/SOURCE/user_data_output_mask.f90
===================================================================
--- palm/trunk/SOURCE/user_data_output_mask.f90	(revision 410)
+++ palm/trunk/SOURCE/user_data_output_mask.f90	(revision 410)
@@ -0,0 +1,72 @@
+ SUBROUTINE user_data_output_mask( av, variable, found, local_pf )
+
+!------------------------------------------------------------------------------!
+! Current revisions:
+! ------------------
+! 
+!
+! Former revisions:
+! -----------------
+! $Id$
+! Initial version
+!
+! Description:
+! ------------
+! Resorts the user-defined output quantity with indices (k,j,i) to a
+! temporary array with indices (i,j,k) for masked data output.
+!------------------------------------------------------------------------------!
+
+    USE control_parameters
+    USE indices
+    USE user
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  variable
+
+    INTEGER ::  av, i, j, k
+
+    LOGICAL ::  found
+
+    REAL, DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2), &
+                    mask_size_l(mid,3)) ::  local_pf
+
+
+    found = .TRUE.
+
+    SELECT CASE ( TRIM( variable ) )
+
+!
+!--    Uncomment and extend the following lines, if necessary.
+!--    The arrays for storing the user defined quantities (here u2 and u2_av)
+!--    have to be declared and defined by the user!
+!--    Sample for user-defined output:
+!       CASE ( 'u2' )
+!          IF ( av == 0 )  THEN
+!            DO  i = 1, mask_size_l(mid,1)
+!               DO  j = 1, mask_size_l(mid,2)
+!                  DO  k = 1, mask_size_l(mid,3)
+!                      local_pf(i,j,k) = u2(mask_k(mid,k), &
+!                                           mask_j(mid,j),mask_i(mid,i))
+!                   ENDDO
+!                ENDDO
+!             ENDDO
+!          ELSE
+!            DO  i = 1, mask_size_l(mid,1)
+!               DO  j = 1, mask_size_l(mid,2)
+!                  DO  k = 1, mask_size_l(mid,3)
+!                      local_pf(i,j,k) = u2_av(mask_k(mid,k), &
+!                                              mask_j(mid,j),mask_i(mid,i))
+!                   ENDDO
+!                ENDDO
+!             ENDDO
+!          ENDIF
+!
+
+       CASE DEFAULT
+          found = .FALSE.
+
+    END SELECT
+
+
+ END SUBROUTINE user_data_output_mask
Index: palm/trunk/SOURCE/user_parin.f90
===================================================================
--- palm/trunk/SOURCE/user_parin.f90	(revision 407)
+++ palm/trunk/SOURCE/user_parin.f90	(revision 410)
@@ -5,4 +5,8 @@
 ! -----------------
 ! 
+!
+! Branch revisions:
+! -----------------
+! masked data output
 !
 ! Former revisions:
@@ -35,5 +39,15 @@
 
 
-    NAMELIST /userpar/  data_output_pr_user, data_output_user, region
+    NAMELIST /userpar/  data_output_pr_user, data_output_user, region,  &
+                        data_output_mask_user_01, data_output_mask_user_02, &
+                        data_output_mask_user_03, data_output_mask_user_04, &
+                        data_output_mask_user_05, data_output_mask_user_06, &
+                        data_output_mask_user_07, data_output_mask_user_08, &
+                        data_output_mask_user_09, data_output_mask_user_10, &
+                        data_output_mask_user_11, data_output_mask_user_12, &
+                        data_output_mask_user_13, data_output_mask_user_14, &
+                        data_output_mask_user_15, data_output_mask_user_16, &
+                        data_output_mask_user_17, data_output_mask_user_18, &
+                        data_output_mask_user_19, data_output_mask_user_20
 
 !
