source: palm/trunk/SOURCE/netcdf_interface_mod.f90 @ 3942

!> @file netcdf_interface_mod.f90
!------------------------------------------------------------------------------!
! This file is part of the PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software
! Foundation, either version 3 of the License, or (at your option) any later
! version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 1997-2019 Leibniz Universitaet Hannover
!------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
! 
! 
! Former revisions:
! -----------------
! $Id: netcdf_interface_mod.f90 3942 2019-04-30 13:08:30Z kanani $
! Add specifier to netcdf_handle_error to simplify identification of attribute
! causing the error
! 
! 3766 2019-02-26 16:23:41Z raasch
! bugfix in im_define_netcdf_grid argument list
! 
! 3745 2019-02-15 18:57:56Z suehring
! Add indoor model
! 
! 3744 2019-02-15 18:38:58Z suehring
! Bugfix: - initialize return values to ensure they are set before returning
!           (routine define_geo_coordinates)
!         - change order of dimensions for some variables
! 
! 3727 2019-02-08 14:52:10Z gronemeier
! make several routines publicly available
! 
! 3701 2019-01-26 18:57:21Z knoop
! Statement added to prevent compiler warning about unused variable
! 
! 3655 2019-01-07 16:51:22Z knoop
! Move the control parameter "salsa" from salsa_mod to control_parameters 
! (M. Kurppa)
! 
! 3582 2018-11-29 19:16:36Z suehring
! dom_dwd_user, Schrempf:
! Remove uv exposure model code, this is now part of biometeorology_mod
! 
! 3529 2018-11-15 21:03:15Z gronemeier
! - set time units
! - add additional global attributes,
! - add additinal variable attributes
! - move definition of UTM and geographic coordinates into subroutine
! - change fill_value
! 
! 3525 2018-11-14 16:06:14Z kanani
! Changes related to clean-up of biometeorology (dom_dwd_user)
! 
! 3485 2018-11-03 17:09:40Z gronemeier
! Write geographic coordinates as global attributes to file.
! 
! 3467 2018-10-30 19:05:21Z suehring
! - Salsa implemented
! - Bugfix convert_utm_to...
! 
! 3464 2018-10-30 18:08:55Z kanani
! - Add variable crs to output files
! - Add long_name to UTM coordinates
! - Add latitude/longitude coordinates. For 3d and xy output, lon and lat are
!   only written if parallel output is used.
! 
! 3459 2018-10-30 15:04:11Z gronemeier
! Adjustment of biometeorology calls
! 
! 3435 2018-10-26 18:25:44Z gronemeier
! Bugfix: corrected order of calls to define_netcdf_grid for masked output
! Add vertical dimensions to masked output in case of terrain-following output
! 
! 3421 2018-10-24 18:39:32Z gronemeier
! Bugfix: move ocean output variables to ocean_mod
! Renamed output variables
! Add UTM coordinates to mask, 3d, xy, xz, yz output
! 
! 3337 2018-10-12 15:17:09Z kanani
! (from branch resler)
! Add biometeorology
! 
! 3294 2018-10-01 02:37:10Z raasch
! changes concerning modularization of ocean option
! 
! 3274 2018-09-24 15:42:55Z knoop
! Modularization of all bulk cloud physics code components
! 
! 3241 2018-09-12 15:02:00Z raasch
! unused variables removed
! 
! 3235 2018-09-07 14:06:15Z sward
! Changed MAS output dimension id_dim_agtnum to be of defined size and no longer
! unlimited. Also changed some MAS output variables to be of type float
! 
! 3198 2018-08-15 09:23:10Z sward
! Redefined MAS limited time dimension to fit usage of multi_agent_system_end
! 
! 3187 2018-07-31 10:32:34Z sward
! Changed agent output to precision NF90_DOUBLE
! 
! 3165 2018-07-24 13:12:42Z sward
! Added agent ID output
! 
! 3159 2018-07-20 11:20:01Z sward
! Added multi agent system
! 
! 3049 2018-05-29 13:52:36Z Giersch
! Error messages revised
! 
! 3045 2018-05-28 07:55:41Z Giersch
! Error messages revised, code adjusted to PALMs coding standards, CASE pt_ext
! pt_new disabled, comment revised
! 
! 3004 2018-04-27 12:33:25Z Giersch
! .NOT. found in if-query added to account for variables found in tcm
! 
! 2964 2018-04-12 16:04:03Z Giersch
! Calculation of fixed number of output time levels for parallel netcdf output
! has been moved completely to check_parameters
! 
! 2932 2018-03-26 09:39:22Z maronga
! Renamed inipar to initialization_parameters.
! 
! 2817 2018-02-19 16:32:21Z knoop
! Preliminary gust module interface implemented
! 
! 2769 2018-01-25 09:22:24Z raasch
! bugfix for calculating number of required output time levels in case of output
! at the beginning of a restart run
! 
! 2766 2018-01-22 17:17:47Z kanani
! Removed preprocessor directive __chem
! 
! 2746 2018-01-15 12:06:04Z suehring
! Move flag plant canopy to modules
! 
! 2718 2018-01-02 08:49:38Z maronga
! Corrected "Former revisions" section
!
! 2696 2017-12-14 17:12:51Z kanani
! Change in file header (GPL part)
! Implementation of uv exposure model (FK)
! Implemented checks for turbulence_closure_mod (TG)
! Implementation of chemistry module (FK)
! Bugfix in setting netcdf grids for LSM variables
! Enable setting of _FillValue attribute in output files (MS)
!
! 2512 2017-10-04 08:26:59Z raasch
! upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny
! no output of ghost layer data any more
! 
! 2302 2017-07-03 14:07:20Z suehring
! Reading of 3D topography using NetCDF data type NC_BYTE
! 
! 2292 2017-06-20 09:51:42Z schwenkel
! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
! includes two more prognostic equations for cloud drop concentration (nc)Â  
! and cloud water content (qc). 
! 
! 2270 2017-06-09 12:18:47Z maronga
! Removed 2 timeseries (shf_eb + qsws_eb). Removed _eb suffixes
! 
! 2265 2017-06-08 16:58:28Z schwenkel
! Unused variables removed.
! 
! 2239 2017-06-01 12:04:51Z suehring
! Bugfix xy-output of land-surface variables
! 
! 2233 2017-05-30 18:08:54Z suehring
!
! 2232 2017-05-30 17:47:52Z suehring
! Adjustments to new topography and surface concept
! 
! Topograpyh height arrays (zu_s_inner, zw_w_inner) are defined locally, output
! only if parallel netcdf. 
!
! Build interface for topography input: 
! - open file in read-only mode
! - read global attributes
! - read variables
!
! Bugfix in xy output (land-surface case)
!
! 2209 2017-04-19 09:34:46Z kanani
! Added support for plant canopy model output
! 
! 2189 2017-03-21 09:29:52Z raasch
! bugfix: rho renamed rho_ocean for the cross section output
!
! 2109 2017-01-10 12:18:08Z raasch
! bugfix: length of character string netcdf_var_name extended to avoid problems
!         which appeared in restart runs due to truncation
!
! 2040 2016-10-26 16:58:09Z gronemeier
! Increased number of possible statistic_regions to 99
! 
! 2037 2016-10-26 11:15:40Z knoop
! Anelastic approximation implemented
! 
! 2031 2016-10-21 15:11:58Z knoop
! renamed variable rho to rho_ocean
! 
! 2011 2016-09-19 17:29:57Z kanani
! Flag urban_surface is now defined in module control_parameters,
! changed prefix for urban surface model output to "usm_",
! introduced control parameter varnamelength for LEN of trimvar.
! 
! 2007 2016-08-24 15:47:17Z kanani
! Added support for new urban surface model (temporary modifications of 
! SELECT CASE ( ) necessary, see variable trimvar),
! increased DIMENSION of do2d_unit, do3d_unit, id_var_do2d, id_var_do3d,
! increased LEN of char_cross_profiles, var_list, var_list_old
! 
! 2000 2016-08-20 18:09:15Z knoop
! Forced header and separation lines into 80 columns
! 
! 1990 2016-08-12 09:54:36Z gronemeier
! Bugfix: variable list was not written for time series output
! 
! 1980 2016-07-29 15:51:57Z suehring
! Bugfix, in order to steer user-defined output, setting flag found explicitly
! to .F.
! 
! 1976 2016-07-27 13:28:04Z maronga
! Removed remaining 2D land surface quantities. Definition of radiation 
! quantities is now done directly in the respective module
! 
! 1972 2016-07-26 07:52:02Z maronga
! Bugfix: wrong units for lsm quantities.
! Definition of grids for land surface quantities is now done directly in the
! respective module.
! 
! 1960 2016-07-12 16:34:24Z suehring
! Additional labels and units for timeseries output of passive scalar-related 
! quantities
! 
! 1957 2016-07-07 10:43:48Z suehring
! flight module added
!
! 1850 2016-04-08 13:29:27Z maronga
! Module renamed
! 
! 
! 1833 2016-04-07 14:23:03Z raasch
! spectrum renamed spectra_mod
!
! 1786 2016-03-08 05:49:27Z raasch
! Bugfix: id_var_time_sp made public
!
! 1783 2016-03-06 18:36:17Z raasch
! netcdf interface has been modularized, former file netcdf renamed to
! netcdf_interface, creation of netcdf-dimensions and -variables moved to
! specific new subroutines create_netcdf_dim and create_netcdf_var,
! compression (deflation) of variables implemented,
! ibmy special cpp directive removed
!
! 1745 2016-02-05 13:06:51Z gronemeier
! Bugfix: recalculating ntdim_3d, ntdim_2d_xy/xz/yz when checking the
!         extensibility of an existing file (only when using parallel NetCDF). 
! 
! 1691 2015-10-26 16:17:44Z maronga
! Added output of radiative heating rates for RRTMG. Corrected output of 
! radiative fluxes
! 
! 1682 2015-10-07 23:56:08Z knoop
! Code annotations made doxygen readable 
! 
! 1596 2015-05-21 09:34:28Z gronemeier
! Bugfix in masked data output. Read 'zu_3d' when trying to extend masked data
! 
! 1551 2015-03-03 14:18:16Z maronga
! Added support for land surface model and radiation model output. In the course
! of this action a new vertical grid zs (soil) was introduced.
! 
! 1353 2014-04-08 15:21:23Z heinze
! REAL constants provided with KIND-attribute
!
! 1322 2014-03-20 16:38:49Z raasch
! Forgotten ONLY-attribute added to USE-statements
!
! 1320 2014-03-20 08:40:49Z raasch
! ONLY-attribute added to USE-statements,
! kind-parameters added to all INTEGER and REAL declaration statements,
! kinds are defined in new module kinds,
! revision history before 2012 removed,
! comment fields (!:) to be used for variable explanations added to
! all variable declaration statements
!
! 1308 2014-03-13 14:58:42Z fricke
! +ntime_count, oldmode
! Adjust NF90_CREATE and NF90_OPEN statement for parallel output
! To increase the performance for parallel output, the following is done:
! - Limit time dimension
! - Values of axis data are only written by PE0
! - No fill is set for all variables
! Check the number of output time levels for restart jobs
!
! 1206 2013-07-18 12:49:16Z witha
! Bugfix: typo in preprocessor directive in subroutine open_write_netcdf_file
!
! 1092 2013-02-02 11:24:22Z raasch
! unused variables removed
!
! 1053 2012-11-13 17:11:03Z hoffmann
! +qr, nr, prr
!
! 1036 2012-10-22 13:43:42Z raasch
! code put under GPL (PALM 3.9)
!
! 1031 2012-10-19 14:35:30Z raasch
! netCDF4 without parallel file support implemented, new routines
! create_netcdf_file and open_write_netcdf_file at end
!
! 992 2012-09-05 15:08:26Z hoffmann
! Removal of the informative messages PA0352 and PA0353.
!  
! 983 2012-08-21 14:17:57Z hoffmann
! Bugfix in cross_profiles.
!
! 964 2012-07-26 09:14:24Z raasch
! rev 951 and 959 reformatted
!
! 959 2012-07-24 13:13:41Z hoffmann
! Bugfix in cross_profiles. It is not allowed to arrange more than 100
! profiles with cross_profiles. 
!
! 951 2012-07-19 14:22:52Z hoffmann
! cross_profiles, profile_rows, profile_columns are written to netCDF header
!
! Revision 1.1Â  2005/05/18 15:37:16Â  raasch
! Initial revision
!
!
! Description:
! ------------
!> In case of extend = .FALSE.:
!> Define all necessary dimensions, axes and variables for the different
!> netCDF datasets. This subroutine is called from check_open after a new
!> dataset is created. It leaves the open netCDF files ready to write.
!>
!> In case of extend = .TRUE.:
!> Find out if dimensions and variables of an existing file match the values
!> of the actual run. If so, get all necessary information (ids, etc.) from
!> this file.
!>
!> Parameter av can assume values 0 (non-averaged data) and 1 (time averaged
!> data)
!>
!> @todo calculation of output time levels for parallel NetCDF still does not
!>Â  Â  Â  Â cover every exception (change of dt_do, end_time in restart)
!> @todo timeseries and profile output still needs to be rewritten to allow
!>Â  Â  Â  Â modularization
!> @todo output 2d UTM coordinates without global arrays
!> @todo output longitude/latitude also with non-parallel output (3d and xy)
!------------------------------------------------------------------------------!
 MODULE netcdf_interface

    USE control_parameters,                                                    &
        ONLY:  biometeorology, fl_max,                                         &
               max_masks, multi_agent_system_end,                              &
               multi_agent_system_start, var_fl_max, varnamelength
    USE kinds
#if defined( __netcdf )
    USE NETCDF
#endif
    USE mas_global_attributes,                                                 &
        ONLY:  dim_size_agtnum

    USE netcdf_data_input_mod,                                                 &
        ONLY: coord_ref_sys, init_model

    PRIVATE

    CHARACTER (LEN=16), DIMENSION(13) ::  agt_var_names =                      &
          (/ 'ag_id           ', 'ag_x            ', 'ag_y            ',       &
             'ag_wind         ', 'ag_temp         ', 'ag_group        ',       &
             'PM10            ', 'PM25            ', 'ag_iPT          ',       &
             'ag_uv           ', 'not_used        ', 'not_used        ',       &
             'not_used        ' /)

    CHARACTER (LEN=16), DIMENSION(13) ::  agt_var_units = &
          (/ 'dim_less        ', 'meters          ', 'meters          ',       &
             'm/s             ', 'K               ', 'dim_less        ',       &
             'tbd             ', 'tbd             ', 'tbd             ',       &
             'tbd             ', 'not_used        ', 'not_used        ',       &
             'not_used        ' /)

    INTEGER(iwp), PARAMETER ::  dopr_norm_num = 7, dopts_num = 29, dots_max = 100

    CHARACTER (LEN=7), DIMENSION(dopr_norm_num) ::  dopr_norm_names =          &
         (/ 'wtheta0', 'ws2    ', 'tsw2   ', 'ws3    ', 'ws2tsw ', 'wstsw2 ',  &
            'z_i    ' /)

    CHARACTER (LEN=7), DIMENSION(dopr_norm_num) ::  dopr_norm_longnames =      &
         (/ 'wtheta0', 'w*2    ', 't*w2   ', 'w*3    ', 'w*2t*w ', 'w*t*w2 ',  &
            'z_i    ' /)

    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_label =                   &
          (/ 'tnpt   ', 'x_     ', 'y_     ', 'z_     ', 'z_abs  ', 'u      ', &
             'v      ', 'w      ', 'u"     ', 'v"     ', 'w"     ', 'npt_up ', &
             'w_up   ', 'w_down ', 'radius ', 'r_min  ', 'r_max  ', 'npt_max', &
             'npt_min', 'x*2    ', 'y*2    ', 'z*2    ', 'u*2    ', 'v*2    ', &
             'w*2    ', 'u"2    ', 'v"2    ', 'w"2    ', 'npt*2  ' /)

    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_unit =                    &
          (/ 'number ', 'm      ', 'm      ', 'm      ', 'm      ', 'm/s    ', &
             'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'number ', &
             'm/s    ', 'm/s    ', 'm      ', 'm      ', 'm      ', 'number ', &
             'number ', 'm2     ', 'm2     ', 'm2     ', 'm2/s2  ', 'm2/s2  ', &
             'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'number2' /)

    INTEGER(iwp) ::  dots_num  = 25  !< number of timeseries defined by default
    INTEGER(iwp) ::  dots_soil = 26  !< starting index for soil-timeseries
    INTEGER(iwp) ::  dots_rad  = 32  !< starting index for radiation-timeseries

    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label =                    &
          (/ 'E            ', 'E*           ', 'dt           ',                &
             'us*          ', 'th*          ', 'umax         ',                &
             'vmax         ', 'wmax         ', 'div_new      ',                &
             'div_old      ', 'zi_wtheta    ', 'zi_theta     ',                &
             'w*           ', 'w"theta"0    ', 'w"theta"     ',                &
             'wtheta       ', 'theta(0)     ', 'theta(z_mo)  ',                &
             'w"u"0        ', 'w"v"0        ', 'w"q"0        ',                &
             'ol           ', 'q*           ', 'w"s"         ',                &
             's*           ', 'ghf          ', 'qsws_liq     ',                &
             'qsws_soil    ', 'qsws_veg     ', 'r_a          ',                &
             'r_s          ',                                                  &
             'rad_net      ', 'rad_lw_in    ', 'rad_lw_out   ',                &
             'rad_sw_in    ', 'rad_sw_out   ', 'rrtm_aldif   ',                &
             'rrtm_aldir   ', 'rrtm_asdif   ', 'rrtm_asdir   ',                &                                               
             ( 'unknown      ', i9 = 1, dots_max-40 ) /)

    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit =                     &
          (/ 'm2/s2        ', 'm2/s2        ', 's            ',                &
             'm/s          ', 'K            ', 'm/s          ',                &
             'm/s          ', 'm/s          ', 's-1          ',                &
             's-1          ', 'm            ', 'm            ',                &
             'm/s          ', 'K m/s        ', 'K m/s        ',                &
             'K m/s        ', 'K            ', 'K            ',                &
             'm2/s2        ', 'm2/s2        ', 'kg m/s       ',                &
             'm            ', 'kg/kg        ', 'kg m/(kg s)  ',                &
             'kg/kg        ', 'W/m2         ', 'W/m2         ',                &
             'W/m2         ', 'W/m2         ', 's/m          ',                &
             's/m          ',                                                  &
             'W/m2         ', 'W/m2         ', 'W/m2         ',                &
             'W/m2         ', 'W/m2         ', '             ',                &
             '             ', '             ', '             ',                &
             ( 'unknown      ', i9 = 1, dots_max-40 ) /)

    CHARACTER (LEN=16) :: heatflux_output_unit     !< unit for heatflux output
    CHARACTER (LEN=16) :: waterflux_output_unit    !< unit for waterflux output
    CHARACTER (LEN=16) :: momentumflux_output_unit !< unit for momentumflux output

    CHARACTER (LEN=9), DIMENSION(300) ::  dopr_unit = 'unknown'

    CHARACTER (LEN=7), DIMENSION(0:1,500) ::  do2d_unit, do3d_unit

!    CHARACTER (LEN=16), DIMENSION(25) ::  prt_var_names = &
!          (/ 'pt_age          ', 'pt_dvrp_size    ', 'pt_origin_x     ', &
!             'pt_origin_y     ', 'pt_origin_z     ', 'pt_radius       ', &
!             'pt_speed_x      ', 'pt_speed_y      ', 'pt_speed_z      ', &
!             'pt_weight_factor', 'pt_x            ', 'pt_y            ', &
!             'pt_z            ', 'pt_color        ', 'pt_group        ', &
!             'pt_tailpoints   ', 'pt_tail_id      ', 'pt_density_ratio', &
!             'pt_exp_arg      ', 'pt_exp_term     ', 'not_used        ', &
!             'not_used        ', 'not_used        ', 'not_used        ', &
!             'not_used        ' /)

!    CHARACTER (LEN=16), DIMENSION(25) ::  prt_var_units = &
!          (/ 'seconds         ', 'meters          ', 'meters          ', &
!             'meters          ', 'meters          ', 'meters          ', &
!             'm/s             ', 'm/s             ', 'm/s             ', &
!             'factor          ', 'meters          ', 'meters          ', &
!             'meters          ', 'none            ', 'none            ', &
!             'none            ', 'none            ', 'ratio           ', &
!             'none            ', 'none            ', 'not_used        ', &
!             'not_used        ', 'not_used        ', 'not_used        ', &
!             'not_used        ' /)

    CHARACTER(LEN=20), DIMENSION(11) ::  netcdf_precision = ' '
    CHARACTER(LEN=40) ::  netcdf_data_format_string

    INTEGER(iwp) ::  id_dim_agtnum, id_dim_time_agt,                           &
                     id_dim_time_fl, id_dim_time_pr,                           &
                     id_dim_time_pts, id_dim_time_sp, id_dim_time_ts,          &
                     id_dim_x_sp, id_dim_y_sp, id_dim_zu_sp, id_dim_zw_sp,     &
                     id_set_agt, id_set_fl, id_set_pr, id_set_prt, id_set_pts, &
                     id_set_sp, id_set_ts, id_var_agtnum, id_var_time_agt,     &
                     id_var_time_fl, id_var_rnoa_agt, id_var_time_pr,          &
                     id_var_time_pts, id_var_time_sp, id_var_time_ts,          &
                     id_var_x_sp, id_var_y_sp, id_var_zu_sp, id_var_zw_sp,     &
                     nc_stat


    INTEGER(iwp), DIMENSION(0:1) ::  id_dim_time_xy, id_dim_time_xz, &
                    id_dim_time_yz, id_dim_time_3d, id_dim_x_xy, id_dim_xu_xy, &
                    id_dim_x_xz, id_dim_xu_xz, id_dim_x_yz, id_dim_xu_yz, &
                    id_dim_x_3d, id_dim_xu_3d, id_dim_y_xy, id_dim_yv_xy, &
                    id_dim_y_xz, id_dim_yv_xz, id_dim_y_yz, id_dim_yv_yz, &
                    id_dim_y_3d, id_dim_yv_3d, id_dim_zs_xy, id_dim_zs_xz, &
                    id_dim_zs_yz, id_dim_zs_3d, id_dim_zu_xy, id_dim_zu1_xy, &
                    id_dim_zu_xz, id_dim_zu_yz, id_dim_zu_3d, id_dim_zw_xy, &
                    id_dim_zw_xz, id_dim_zw_yz, id_dim_zw_3d, id_set_xy, &
                    id_set_xz, id_set_yz, id_set_3d, id_var_ind_x_yz, &
                    id_var_ind_y_xz, id_var_ind_z_xy, id_var_time_xy, &
                    id_var_time_xz, id_var_time_yz, id_var_time_3d, id_var_x_xy, &
                    id_var_xu_xy, id_var_x_xz, id_var_xu_xz, id_var_x_yz, &
                    id_var_xu_yz, id_var_x_3d, id_var_xu_3d, id_var_y_xy, &
                    id_var_yv_xy, id_var_y_xz, id_var_yv_xz, id_var_y_yz, &
                    id_var_yv_yz, id_var_y_3d, id_var_yv_3d, id_var_zs_xy, &
                    id_var_zs_xz, id_var_zs_yz, id_var_zs_3d, id_var_zusi_xy, &
                    id_var_zusi_3d, id_var_zu_xy, id_var_zu1_xy, id_var_zu_xz, &
                    id_var_zu_yz, id_var_zu_3d, id_var_zwwi_xy, id_var_zwwi_3d, &
                    id_var_zw_xy, id_var_zw_xz, id_var_zw_yz, id_var_zw_3d

    INTEGER(iwp), DIMENSION(0:2,0:1) ::  id_var_eutm_3d, id_var_nutm_3d, &
                                         id_var_eutm_xy, id_var_nutm_xy, &
                                         id_var_eutm_xz, id_var_nutm_xz, &
                                         id_var_eutm_yz, id_var_nutm_yz

    INTEGER(iwp), DIMENSION(0:2,0:1) ::  id_var_lat_3d, id_var_lon_3d, &
                                         id_var_lat_xy, id_var_lon_xy, &
                                         id_var_lat_xz, id_var_lon_xz, &
                                         id_var_lat_yz, id_var_lon_yz

    INTEGER ::  netcdf_data_format = 2  !< NetCDF3 64bit offset format
    INTEGER ::  netcdf_deflate = 0      !< NetCDF compression, default: no
                                        !< compression

    INTEGER(iwp)                 ::  dofl_time_count
    INTEGER(iwp), DIMENSION(10)  ::  id_var_dospx, id_var_dospy
    INTEGER(iwp), DIMENSION(20)  ::  id_var_agt
!    INTEGER(iwp), DIMENSION(20)  ::  id_var_prt
    INTEGER(iwp), DIMENSION(11)  ::  nc_precision
    INTEGER(iwp), DIMENSION(dopr_norm_num) ::  id_var_norm_dopr
    
    INTEGER(iwp), DIMENSION(fl_max) ::  id_dim_x_fl, id_dim_y_fl, id_dim_z_fl
    INTEGER(iwp), DIMENSION(fl_max) ::  id_var_x_fl, id_var_y_fl, id_var_z_fl
    
    CHARACTER (LEN=20), DIMENSION(fl_max*var_fl_max) :: dofl_label
    CHARACTER (LEN=20), DIMENSION(fl_max*var_fl_max) :: dofl_unit 
    CHARACTER (LEN=20), DIMENSION(fl_max) :: dofl_dim_label_x
    CHARACTER (LEN=20), DIMENSION(fl_max) :: dofl_dim_label_y
    CHARACTER (LEN=20), DIMENSION(fl_max) :: dofl_dim_label_z

    INTEGER(iwp), DIMENSION(fl_max*var_fl_max) :: id_var_dofl    

    INTEGER(iwp), DIMENSION(dopts_num,0:10) ::  id_var_dopts
    INTEGER(iwp), DIMENSION(0:1,500)        ::  id_var_do2d, id_var_do3d
    INTEGER(iwp), DIMENSION(100,0:99)       ::  id_dim_z_pr, id_var_dopr, &
                                                id_var_z_pr
    INTEGER(iwp), DIMENSION(dots_max,0:99)  ::  id_var_dots

!
!-- Masked output
    CHARACTER (LEN=7), DIMENSION(max_masks,0:1,100) ::  domask_unit

    LOGICAL ::  output_for_t0 = .FALSE.

    INTEGER(iwp), 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_zs_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_zs_mask, &
                   id_var_zu_mask, id_var_zw_mask, &
                   id_var_zusi_mask, id_var_zwwi_mask

    INTEGER(iwp), DIMENSION(0:2,1:max_masks,0:1) ::  id_var_eutm_mask, &
                                                     id_var_nutm_mask

    INTEGER(iwp), DIMENSION(0:2,1:max_masks,0:1) ::  id_var_lat_mask, &
                                                     id_var_lon_mask

    INTEGER(iwp), DIMENSION(1:max_masks,0:1,100) ::  id_var_domask

    REAL(wp) ::  fill_value = -9999.0_wp    !< value for the _FillValue attribute


    PUBLIC  dofl_dim_label_x, dofl_dim_label_y, dofl_dim_label_z, dofl_label,  &
            dofl_time_count, dofl_unit, domask_unit, dopr_unit, dopts_num,     &
            dots_label, dots_max, dots_num, dots_rad, dots_soil, dots_unit,    &
            do2d_unit, do3d_unit, fill_value, id_set_agt, id_set_fl,           &
            id_set_mask, id_set_pr, id_set_prt, id_set_pts, id_set_sp,         &
            id_set_ts, id_set_xy, id_set_xz, id_set_yz, id_set_3d, id_var_agt, &
            id_var_domask, id_var_dofl, id_var_dopr, id_var_dopts,             &
            id_var_dospx, id_var_dospy, id_var_dots, id_var_do2d, id_var_do3d, &
            id_var_norm_dopr, id_var_time_agt, id_var_time_fl,                 &
            id_var_time_mask, id_var_time_pr, id_var_rnoa_agt, id_var_time_pts,&
            id_var_time_sp, id_var_time_ts,                                    &
            id_var_time_xy, id_var_time_xz, id_var_time_yz, id_var_time_3d,    &
            id_var_x_fl, id_var_y_fl, id_var_z_fl,  nc_stat,                   &
            netcdf_data_format, netcdf_data_format_string, netcdf_deflate,     &
            netcdf_precision, output_for_t0, heatflux_output_unit,             &
            waterflux_output_unit, momentumflux_output_unit

    SAVE

    INTERFACE netcdf_create_dim
       MODULE PROCEDURE netcdf_create_dim
    END INTERFACE netcdf_create_dim

    INTERFACE netcdf_create_file
       MODULE PROCEDURE netcdf_create_file
    END INTERFACE netcdf_create_file

    INTERFACE netcdf_create_global_atts
       MODULE PROCEDURE netcdf_create_global_atts
    END INTERFACE netcdf_create_global_atts

    INTERFACE netcdf_create_var
       MODULE PROCEDURE netcdf_create_var
    END INTERFACE netcdf_create_var

    INTERFACE netcdf_create_att
       MODULE PROCEDURE netcdf_create_att_string
    END INTERFACE netcdf_create_att

    INTERFACE netcdf_define_header
       MODULE PROCEDURE netcdf_define_header
    END INTERFACE netcdf_define_header

    INTERFACE netcdf_handle_error
       MODULE PROCEDURE netcdf_handle_error
    END INTERFACE netcdf_handle_error

    INTERFACE netcdf_open_write_file
       MODULE PROCEDURE netcdf_open_write_file
    END INTERFACE netcdf_open_write_file

    PUBLIC netcdf_create_att, netcdf_create_dim, netcdf_create_file,           &
           netcdf_create_global_atts, netcdf_create_var, netcdf_define_header, &
           netcdf_handle_error, netcdf_open_write_file

 CONTAINS

 SUBROUTINE netcdf_define_header( callmode, extend, av )
 
#if defined( __netcdf )

    USE arrays_3d,                                                             &
        ONLY:  zu, zw

    USE biometeorology_mod,                                                    &
        ONLY:  bio_define_netcdf_grid

    USE chemistry_model_mod,                                                   &
        ONLY:  chem_define_netcdf_grid 

    USE basic_constants_and_equations_mod,                                     &
        ONLY:  pi

    USE control_parameters,                                                    &
        ONLY:  agent_time_unlimited, air_chemistry, averaging_interval,        &
               averaging_interval_pr, data_output_pr, domask, dopr_n,          &
               dopr_time_count, dopts_time_count, dots_time_count,             &
               do2d, do2d_at_begin, do2d_xz_time_count, do3d, do3d_at_begin,   &
               do2d_yz_time_count, dt_data_output_av, dt_do2d_xy, dt_do2d_xz,  &
               dt_do2d_yz, dt_do3d, dt_write_agent_data, mask_size,            &
               do2d_xy_time_count, do3d_time_count, domask_time_count,         &
               end_time, indoor_model, land_surface,                           &
               mask_size_l, mask_i, mask_i_global, mask_j, mask_j_global,      &
               mask_k_global, mask_surface,                                    &
               message_string, mid, ntdim_2d_xy, ntdim_2d_xz,                  &
               ntdim_2d_yz, ntdim_3d, nz_do3d, ocean_mode, plant_canopy,       &
               run_description_header, salsa, section, simulated_time,         &
               simulated_time_at_begin, skip_time_data_output_av,              &
               skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,        &
               skip_time_do3d, topography, num_leg, num_var_fl,                &
               urban_surface

    USE grid_variables,                                                        &
        ONLY:  dx, dy, zu_s_inner, zw_w_inner

    USE gust_mod,                                                              &
        ONLY: gust_define_netcdf_grid, gust_module_enabled

    USE indices,                                                               &
        ONLY:  nx, nxl, nxr, ny, nys, nyn, nz ,nzb, nzt

    USE kinds

    USE indoor_model_mod,                                                      &
        ONLY: im_define_netcdf_grid
    
    USE land_surface_model_mod,                                                &
        ONLY: lsm_define_netcdf_grid, nzb_soil, nzt_soil, nzs, zs

    USE ocean_mod,                                                             &
        ONLY:  ocean_define_netcdf_grid

    USE pegrid

    USE particle_attributes,                                                   &
        ONLY:  number_of_particle_groups

    USE plant_canopy_model_mod,                                                &
        ONLY:  pcm_define_netcdf_grid

    USE profil_parameter,                                                      &
        ONLY:  crmax, cross_profiles, dopr_index, profile_columns, profile_rows

    USE radiation_model_mod,                                                   &
        ONLY: radiation, radiation_define_netcdf_grid      
     
    USE salsa_mod,                                                             &
        ONLY:  salsa_define_netcdf_grid            

    USE spectra_mod,                                                           &
        ONLY:  averaging_interval_sp, comp_spectra_level, data_output_sp, dosp_time_count, spectra_direction

    USE statistics,                                                            &
        ONLY:  hom, statistic_regions

    USE turbulence_closure_mod,                                                &
        ONLY:  tcm_define_netcdf_grid

    USE urban_surface_mod,                                                     &
        ONLY:  usm_define_netcdf_grid

    USE user,                                                                  &
        ONLY:  user_module_enabled, user_define_netcdf_grid



    IMPLICIT NONE

    CHARACTER (LEN=3)              ::  suffix                !<
    CHARACTER (LEN=2), INTENT (IN) ::  callmode              !<
    CHARACTER (LEN=4)              ::  grid_x                !<
    CHARACTER (LEN=4)              ::  grid_y                !<
    CHARACTER (LEN=4)              ::  grid_z                !<
    CHARACTER (LEN=6)              ::  mode                  !<
    CHARACTER (LEN=10)             ::  precision             !<
    CHARACTER (LEN=10)             ::  var                   !<
    CHARACTER (LEN=20)             ::  netcdf_var_name       !<
    CHARACTER (LEN=varnamelength)  ::  trimvar               !< TRIM of output-variable string
    CHARACTER (LEN=80)             ::  time_average_text     !<
    CHARACTER (LEN=4000)           ::  char_cross_profiles   !<
    CHARACTER (LEN=4000)           ::  var_list              !<
    CHARACTER (LEN=4000)           ::  var_list_old          !<

    CHARACTER (LEN=100), DIMENSION(1:crmax) ::  cross_profiles_adj   !<
    CHARACTER (LEN=100), DIMENSION(1:crmax) ::  cross_profiles_char  !<

    INTEGER(iwp) ::  av                                      !<
    INTEGER(iwp) ::  cross_profiles_count                    !<
    INTEGER(iwp) ::  cross_profiles_maxi                     !<
    INTEGER(iwp) ::  delim                                   !<
    INTEGER(iwp) ::  delim_old                               !<
    INTEGER(iwp) ::  file_id                                 !<
    INTEGER(iwp) ::  i                                       !<
    INTEGER(iwp) ::  id_last                                 !<
    INTEGER(iwp) ::  id_x                                    !<
    INTEGER(iwp) ::  id_y                                    !<
    INTEGER(iwp) ::  id_z                                    !<
    INTEGER(iwp) ::  j                                       !<
    INTEGER(iwp) ::  k                                       !<
    INTEGER(iwp) ::  kk                                      !<
    INTEGER(iwp) ::  ns                                      !<
    INTEGER(iwp) ::  ns_do                                   !< actual value of ns for soil model data
    INTEGER(iwp) ::  ns_old                                  !<
    INTEGER(iwp) ::  ntime_count                             !< number of time levels found in file
    INTEGER(iwp) ::  nz_old                                  !<
    INTEGER(iwp) ::  l                                       !<

    INTEGER(iwp), SAVE ::  oldmode                           !<

    INTEGER(iwp), DIMENSION(1) ::  id_dim_time_old           !<
    INTEGER(iwp), DIMENSION(1) ::  id_dim_x_yz_old           !<
    INTEGER(iwp), DIMENSION(1) ::  id_dim_y_xz_old           !<
    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_sp_old          !<
    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_xy_old          !<
    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_3d_old          !<
    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_mask_old        !<


    INTEGER(iwp), DIMENSION(1:crmax) ::  cross_profiles_numb !<

    LOGICAL ::  found                                        !<

    LOGICAL, INTENT (INOUT) ::  extend                       !<

    LOGICAL, SAVE ::  init_netcdf = .FALSE.                  !<

    REAL(wp) ::  cos_ra                                      !< cosine of rotation_angle
    REAL(wp) ::  eutm                                        !< easting (UTM)
    REAL(wp) ::  nutm                                        !< northing (UTM)
    REAL(wp) ::  shift_x                                     !< shift of x coordinate
    REAL(wp) ::  shift_y                                     !< shift of y coordinate
    REAL(wp) ::  sin_ra                                      !< sine of rotation_angle

    REAL(wp), DIMENSION(1) ::  last_time_coordinate          !< last time value in file
    REAL(wp), DIMENSION(8) ::  crs_list                      !< list of coord_ref_sys values

    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  netcdf_data    !<
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  netcdf_data_2d !<
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lat            !< latitude
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lon            !< longitude


!
!-- Initializing actions
    IF ( .NOT. init_netcdf )  THEN
!
!--Â  Â  Check and set accuracy for netCDF output. First set default value
       nc_precision = NF90_REAL4

       i = 1
       DO  WHILE ( netcdf_precision(i) /= ' ' )
          j = INDEX( netcdf_precision(i), '_' )
          IF ( j == 0 )  THEN
             WRITE ( message_string, * ) 'netcdf_precision must contain a ', &
                                         '"_"netcdf_precision(', i, ')="',   &
                                         TRIM( netcdf_precision(i) ),'"'
             CALL message( 'netcdf_define_header', 'PA0241', 2, 2, 0, 6, 0 )
          ENDIF

          var       = netcdf_precision(i)(1:j-1)
          precision = netcdf_precision(i)(j+1:)

          IF ( precision == 'NF90_REAL4' )  THEN
             j = NF90_REAL4
          ELSEIF ( precision == 'NF90_REAL8' )  THEN
             j = NF90_REAL8
          ELSE
             WRITE ( message_string, * ) 'illegal netcdf precision: ',  &
                                         'netcdf_precision(', i, ')="', &
                                         TRIM( netcdf_precision(i) ),'"'
             CALL message( 'netcdf_define_header', 'PA0242', 1, 2, 0, 6, 0 )
          ENDIF

          SELECT CASE ( var )
             CASE ( 'xy' )
                nc_precision(1) = j
             CASE ( 'xz' )
                nc_precision(2) = j
             CASE ( 'yz' )
                nc_precision(3) = j
             CASE ( '2d' )
                nc_precision(1:3) = j
             CASE ( '3d' )
                nc_precision(4) = j
             CASE ( 'pr' )
                nc_precision(5) = j
             CASE ( 'ts' )
                nc_precision(6) = j
             CASE ( 'sp' )
                nc_precision(7) = j
             CASE ( 'prt' )
                nc_precision(8) = j
             CASE ( 'masks' )
                nc_precision(11) = j
             CASE ( 'fl' )
                nc_precision(9) = j
             CASE ( 'all' )
                nc_precision    = j

             CASE DEFAULT
                WRITE ( message_string, * ) 'unknown variable in ' //          &
                                  'initialization_parameters ',                & 
                                  'assignment: netcdf_precision(', i, ')="',   &
                                            TRIM( netcdf_precision(i) ),'"'
                CALL message( 'netcdf_define_header', 'PA0243', 1, 2, 0, 6, 0 )

          END SELECT

          i = i + 1
          IF ( i > 50 )  EXIT
       ENDDO

!
!--Â  Â  Check for allowed parameter range
       IF ( netcdf_deflate < 0  .OR.  netcdf_deflate > 9 )  THEN
          WRITE ( message_string, '(A,I3,A)' ) 'netcdf_deflate out of ' //     &
                                      'range & given value: ', netcdf_deflate, &
                                      ', allowed range: 0-9'
          CALL message( 'netcdf_define_header', 'PA0355', 2, 2, 0, 6, 0 )
       ENDIF
!
!--Â  Â  Data compression does not work with parallel NetCDF/HDF5
       IF ( netcdf_deflate > 0  .AND.  netcdf_data_format /= 3 )  THEN
          message_string = 'netcdf_deflate reset to 0'
          CALL message( 'netcdf_define_header', 'PA0356', 0, 1, 0, 6, 0 )

          netcdf_deflate = 0
       ENDIF

       init_netcdf = .TRUE.

    ENDIF

!
!-- Convert coord_ref_sys into vector (used for lat/lon calculation)
    crs_list = (/ coord_ref_sys%semi_major_axis,                  &
                  coord_ref_sys%inverse_flattening,               &
                  coord_ref_sys%longitude_of_prime_meridian,      &
                  coord_ref_sys%longitude_of_central_meridian,    &
                  coord_ref_sys%scale_factor_at_central_meridian, &
                  coord_ref_sys%latitude_of_projection_origin,    &
                  coord_ref_sys%false_easting,                    &
                  coord_ref_sys%false_northing /)

!
!-- Determine the mode to be processed
    IF ( extend )  THEN
       mode = callmode // '_ext'
    ELSE
       mode = callmode // '_new'
    ENDIF

!
!-- Select the mode to be processed. Possibilities are 3d, ma (mask), xy, xz, 
!-- yz, pr (profiles), ps (particle timeseries), fl (flight data), ts 
!-- (timeseries) or sp (spectra)
    SELECT CASE ( mode )

       CASE ( 'ma_new' )

!
!--Â  Â  Â  Â decompose actual parameter file_id (=formal parameter av) into 
!--Â  Â  Â  Â mid and av
          file_id = av
          IF ( file_id <= 200+max_masks )  THEN
             mid = file_id - 200
             av = 0
          ELSE
             mid = file_id - (200+max_masks)
             av = 1
          ENDIF

!
!--Â  Â  Â  Â Define some global attributes of the dataset
          IF ( av == 0 )  THEN
             CALL netcdf_create_global_atts( id_set_mask(mid,av), 'podsmasked', TRIM( run_description_header ), 464 )
             time_average_text = ' '
          ELSE
             CALL netcdf_create_global_atts( id_set_mask(mid,av), 'podsmasked', TRIM( run_description_header ), 464 )
             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 netcdf_handle_error( 'netcdf_define_header', 466 )
          ENDIF

!
!--Â  Â  Â  Â Define time coordinate for volume data (unlimited dimension)
          CALL netcdf_create_dim( id_set_mask(mid,av), 'time', NF90_UNLIMITED, &
                                  id_dim_time_mask(mid,av), 467 )
          CALL netcdf_create_var( id_set_mask(mid,av),                         &
                                  (/ id_dim_time_mask(mid,av) /), 'time',      &
                                  NF90_DOUBLE, id_var_time_mask(mid,av),       &
                                 'seconds since '//TRIM(init_model%origin_time), 'time', 468, 469, 000 )
          CALL netcdf_create_att( id_set_mask(mid,av), id_var_time_mask(mid,av), 'standard_name', 'time', 000)
          CALL netcdf_create_att( id_set_mask(mid,av), id_var_time_mask(mid,av), 'axis', 'T', 000)

!
!--Â  Â  Â  Â Define spatial dimensions and coordinates:
          IF ( mask_surface(mid) )  THEN
!
!--Â  Â  Â  Â  Â  In case of terrain-following output, the vertical dimensions are
!--Â  Â  Â  Â  Â  indices, not meters
             CALL netcdf_create_dim( id_set_mask(mid,av), 'ku_above_surf',     &
                                     mask_size(mid,3), id_dim_zu_mask(mid,av), &
                                     470 )
             CALL netcdf_create_var( id_set_mask(mid,av),                      &
                                     (/ id_dim_zu_mask(mid,av) /),             &
                                     'ku_above_surf',                          &
                                     NF90_DOUBLE, id_var_zu_mask(mid,av),      &
                                     '1', 'grid point above terrain',          &
                                     471, 472, 000 )
             CALL netcdf_create_dim( id_set_mask(mid,av), 'kw_above_surf',     &
                                     mask_size(mid,3), id_dim_zw_mask(mid,av), &
                                     473 )
             CALL netcdf_create_var( id_set_mask(mid,av),                      &
                                     (/ id_dim_zw_mask(mid,av) /),             &
                                     'kw_above_surf',                          &
                                     NF90_DOUBLE, id_var_zw_mask(mid,av),      &
                                    '1', 'grid point above terrain',           &
                                    474, 475, 000 )
          ELSE
!
!--Â  Â  Â  Â  Â  Define vertical coordinate grid (zu grid)
             CALL netcdf_create_dim( id_set_mask(mid,av), 'zu_3d',             &
                                     mask_size(mid,3), id_dim_zu_mask(mid,av), &
                                     470 )
             CALL netcdf_create_var( id_set_mask(mid,av),                      &
                                     (/ id_dim_zu_mask(mid,av) /), 'zu_3d',    &
                                     NF90_DOUBLE, id_var_zu_mask(mid,av),      &
                                     'meters', '', 471, 472, 000 )
!
!--Â  Â  Â  Â  Â  Define vertical coordinate grid (zw grid)
             CALL netcdf_create_dim( id_set_mask(mid,av), 'zw_3d',             &
                                     mask_size(mid,3), id_dim_zw_mask(mid,av), &
                                     473 )
             CALL netcdf_create_var( id_set_mask(mid,av),                      &
                                     (/ id_dim_zw_mask(mid,av) /), 'zw_3d',    &
                                     NF90_DOUBLE, id_var_zw_mask(mid,av),      &
                                    'meters', '', 474, 475, 000 )
          ENDIF
!
!--Â  Â  Â  Â Define x-axis (for scalar position)
          CALL netcdf_create_dim( id_set_mask(mid,av), 'x', mask_size(mid,1),  &
                                  id_dim_x_mask(mid,av), 476 )
          CALL netcdf_create_var( id_set_mask(mid,av),                         &
                                  (/ id_dim_x_mask(mid,av) /), 'x',            &
                                  NF90_DOUBLE, id_var_x_mask(mid,av),          &
                                  'meters', '', 477, 478, 000 )
!
!--Â  Â  Â  Â Define x-axis (for u position)
          CALL netcdf_create_dim( id_set_mask(mid,av), 'xu', mask_size(mid,1), &
                                  id_dim_xu_mask(mid,av), 479 )
          CALL netcdf_create_var( id_set_mask(mid,av),                         &
                                  (/ id_dim_xu_mask(mid,av) /), 'xu',          &
                                  NF90_DOUBLE, id_var_xu_mask(mid,av),         &
                                  'meters', '', 480, 481, 000 )
!
!--Â  Â  Â  Â Define y-axis (for scalar position)
          CALL netcdf_create_dim( id_set_mask(mid,av), 'y', mask_size(mid,2),  &
                                  id_dim_y_mask(mid,av), 482 )
          CALL netcdf_create_var( id_set_mask(mid,av),                         &
                                  (/ id_dim_y_mask(mid,av) /), 'y',            &
                                  NF90_DOUBLE, id_var_y_mask(mid,av),          &
                                  'meters', '', 483, 484, 000 )
!
!--Â  Â  Â  Â Define y-axis (for v position)
          CALL netcdf_create_dim( id_set_mask(mid,av), 'yv', mask_size(mid,2), &
                                  id_dim_yv_mask(mid,av), 485 )
          CALL netcdf_create_var( id_set_mask(mid,av),                         &
                                  (/ id_dim_yv_mask(mid,av) /),                &
                                  'yv', NF90_DOUBLE, id_var_yv_mask(mid,av),   &
                                  'meters', '', 486, 487, 000 )
!
!--Â  Â  Â  Â Define UTM and geographic coordinates
          CALL define_geo_coordinates( id_set_mask(mid,av),               &
                  (/ id_dim_x_mask(mid,av), id_dim_xu_mask(mid,av) /),    &
                  (/ id_dim_y_mask(mid,av), id_dim_yv_mask(mid,av) /),    &
                  id_var_eutm_mask(:,mid,av), id_var_nutm_mask(:,mid,av), &
                  id_var_lat_mask(:,mid,av), id_var_lon_mask(:,mid,av)    )
!
!--Â  Â  Â  Â Define coordinate-reference system
          CALL netcdf_create_crs( id_set_mask(mid,av), 000 )
!
!--Â  Â  Â  Â In case of non-flat topography define 2d-arrays containing the height
!--Â  Â  Â  Â information. Only for parallel netcdf output.
          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
               netcdf_data_format > 4 )  THEN
!
!--Â  Â  Â  Â  Â  Define zusi = zu(nzb_s_inner)
             CALL netcdf_create_var( id_set_mask(mid,av),                      &
                                     (/ id_dim_x_mask(mid,av),                 &
                                        id_dim_y_mask(mid,av) /), 'zusi',      &
                                     NF90_DOUBLE, id_var_zusi_mask(mid,av),    &
                                     'meters', 'zu(nzb_s_inner)', 488, 489,    &
                                     490 )
!             
!--Â  Â  Â  Â  Â  Define zwwi = zw(nzb_w_inner)
             CALL netcdf_create_var( id_set_mask(mid,av),                      &
                                     (/ id_dim_x_mask(mid,av),                 &
                                        id_dim_y_mask(mid,av) /), 'zwwi',      &
                                     NF90_DOUBLE, id_var_zwwi_mask(mid,av),    &
                                     'meters', 'zw(nzb_w_inner)', 491, 492,    &
                                     493 )
          ENDIF             
 
          IF ( land_surface )  THEN
!
!--Â  Â  Â  Â  Â  Define vertical coordinate grid (zw grid)
             CALL netcdf_create_dim( id_set_mask(mid,av), 'zs_3d',             &
                                     mask_size(mid,3), id_dim_zs_mask(mid,av), &
                                     536 )
             CALL netcdf_create_var( id_set_mask(mid,av),                      &
                                     (/ id_dim_zs_mask(mid,av) /), 'zs_3d',    &
                                     NF90_DOUBLE, id_var_zs_mask(mid,av),      &
                                     'meters', '', 537, 555, 000 )
          ENDIF

!
!--Â  Â  Â  Â Define the variables
          var_list = ';'
          i = 1

          DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
!
!--Â  Â  Â  Â  Â  Temporary solution to account for data output within the new urban 
!--Â  Â  Â  Â  Â  surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
             trimvar = TRIM( domask(mid,av,i) )
             IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
                trimvar = 'usm_output'
             ENDIF
!
!--Â  Â  Â  Â  Â  Check for the grid
             found = .FALSE.
             SELECT CASE ( trimvar )
!
!--Â  Â  Â  Â  Â  Â  Â Most variables are defined on the scalar grid
                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',                &
                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
                       's', 'theta', 'thetal', 'thetav' )

                   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'

!             
!--Â  Â  Â  Â  Â  Â  Â Block of urban surface model outputs
                CASE ( 'usm_output' )

                   CALL usm_define_netcdf_grid( domask( mid,av,i), found,      &
                                                        grid_x, grid_y, grid_z )

                CASE DEFAULT
!
!--                Check for quantities defined in other modules                                                        
                   CALL tcm_define_netcdf_grid( domask( mid,av,i), found,      &
                                                        grid_x, grid_y, grid_z )

                   IF ( .NOT. found  .AND.  air_chemistry )  THEN
                      CALL chem_define_netcdf_grid( domask(mid,av,i), found,   &
                                                    grid_x, grid_y, grid_z )
                   ENDIF

                   IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
                      CALL gust_define_netcdf_grid( domask(mid,av,i), found,   &
                                                    grid_x, grid_y, grid_z )
                   ENDIF

                   IF ( land_surface )  THEN
                      CALL lsm_define_netcdf_grid( domask(mid,av,i), found,    &
                                                   grid_x, grid_y, grid_z )
                   ENDIF

                   IF ( .NOT. found  .AND.  ocean_mode )  THEN
                      CALL ocean_define_netcdf_grid( domask(mid,av,i), found,  &
                                                     grid_x, grid_y, grid_z )
                   ENDIF

                   IF ( .NOT. found  .AND.  plant_canopy )  THEN
                      CALL pcm_define_netcdf_grid( domask(mid,av,i), found,    &
                                                   grid_x, grid_y, grid_z )
                   ENDIF

                   IF ( .NOT. found  .AND.  radiation )  THEN
                      CALL radiation_define_netcdf_grid( domask(mid,av,i),     &
                                                         found, grid_x, grid_y,&
                                                         grid_z )
                   ENDIF
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Check for SALSA quantities
                   IF ( .NOT. found  .AND.  salsa )  THEN
                      CALL salsa_define_netcdf_grid( domask(mid,av,i), found,  &
                                                     grid_x, grid_y, grid_z )
                   ENDIF                                   
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Now check for user-defined quantities
                   IF ( .NOT. found  .AND.  user_module_enabled )  THEN
                      CALL user_define_netcdf_grid( domask(mid,av,i), found,   &
                                                    grid_x, grid_y, grid_z )
                   ENDIF
                                                 
                   IF ( .NOT. found )  THEN
                      WRITE ( message_string, * ) 'no grid defined for',       &
                           ' variable ', TRIM( domask(mid,av,i) )
                      CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, &
                                    6, 0 )
                   ENDIF

             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)
             ELSEIF ( grid_z == "zs" )  THEN
                id_z = id_dim_zs_mask(mid,av)
             ENDIF

!
!--Â  Â  Â  Â  Â  Define the grid
             CALL netcdf_create_var( id_set_mask(mid,av), (/ id_x, id_y, id_z, &
                                     id_dim_time_mask(mid,av) /),              &
                                     domask(mid,av,i), nc_precision(11),       &
                                     id_var_domask(mid,av,i),                  &
                                     TRIM( domask_unit(mid,av,i) ),            &
                                     domask(mid,av,i), 494, 495, 496, .TRUE. )

             var_list = TRIM( var_list ) // TRIM( domask(mid,av,i) ) // ';'

             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 netcdf_handle_error( 'netcdf_define_header', 497 )

!
!--Â  Â  Â  Â Leave netCDF define mode
          nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 498 )

!
!--Â  Â  Â  Â 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_wp ) * 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 netcdf_handle_error( 'netcdf_define_header', 499 )

          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 netcdf_handle_error( 'netcdf_define_header', 500 )

          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_wp ) * 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 netcdf_handle_error( 'netcdf_define_header', 501 )

          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 netcdf_handle_error( 'netcdf_define_header', 502 )

          DEALLOCATE( netcdf_data )

!
!--Â  Â  Â  Â Write UTM coordinates
          IF ( init_model%rotation_angle == 0.0_wp )  THEN
!
!--Â  Â  Â  Â  Â  1D in case of no rotation
             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
!
!--Â  Â  Â  Â  Â  x coordinates
             ALLOCATE( netcdf_data(mask_size(mid,1)) )
             DO  k = 0, 2
!            
!--Â  Â  Â  Â  Â  Â  Â Scalar grid points
                IF ( k == 0 )  THEN
                   shift_x = 0.5
!            
!--Â  Â  Â  Â  Â  Â  Â u grid points
                ELSEIF ( k == 1 )  THEN
                   shift_x = 0.0
!            
!--Â  Â  Â  Â  Â  Â  Â v grid points
                ELSEIF ( k == 2 )  THEN
                   shift_x = 0.5
                ENDIF

                netcdf_data = init_model%origin_x + cos_ra                     &
                       * ( mask_i_global(mid,:mask_size(mid,1)) + shift_x ) * dx

                nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
                                        id_var_eutm_mask(k,mid,av), &
                                        netcdf_data, start = (/ 1 /), &
                                        count = (/ mask_size(mid,1) /) )
                CALL netcdf_handle_error( 'netcdf_define_header', 555 )

             ENDDO
             DEALLOCATE( netcdf_data )
!
!--Â  Â  Â  Â  Â  y coordinates
             ALLOCATE( netcdf_data(mask_size(mid,2)) )
             DO  k = 0, 2
!
!--Â  Â  Â  Â  Â  Â  Â Scalar grid points
                IF ( k == 0 )  THEN
                   shift_y = 0.5
!
!--Â  Â  Â  Â  Â  Â  Â u grid points
                ELSEIF ( k == 1 )  THEN
                   shift_y = 0.5
!
!--Â  Â  Â  Â  Â  Â  Â v grid points
                ELSEIF ( k == 2 )  THEN
                   shift_y = 0.0
                ENDIF

                netcdf_data = init_model%origin_y + cos_ra                     &
                       * ( mask_j_global(mid,:mask_size(mid,2)) + shift_y ) * dy

                nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
                                        id_var_nutm_mask(k,mid,av), &
                                        netcdf_data, start = (/ 1 /), &
                                        count = (/ mask_size(mid,2) /) )
                CALL netcdf_handle_error( 'netcdf_define_header', 556 )

             ENDDO
             DEALLOCATE( netcdf_data )

          ELSE
!
!--Â  Â  Â  Â  Â  2D in case of rotation
             ALLOCATE( netcdf_data_2d(mask_size(mid,1),mask_size(mid,2)) )
             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
             sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )

             DO  k = 0, 2
!            
!--Â  Â  Â  Â  Â  Â  Scalar grid points
               IF ( k == 0 )  THEN
                  shift_x = 0.5 ; shift_y = 0.5
!            
!--Â  Â  Â  Â  Â  Â  u grid points
               ELSEIF ( k == 1 )  THEN
                  shift_x = 0.0 ; shift_y = 0.5
!            
!--Â  Â  Â  Â  Â  Â  v grid points
               ELSEIF ( k == 2 )  THEN
                  shift_x = 0.5 ; shift_y = 0.0
               ENDIF

               DO  j = 1, mask_size(mid,2)
                  DO  i = 1, mask_size(mid,1)
                     netcdf_data_2d(i,j) = init_model%origin_x                 &
                           + cos_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
                           + sin_ra * ( mask_j_global(mid,j) + shift_y ) * dy
                  ENDDO
               ENDDO

               nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
                                       id_var_eutm_mask(k,mid,av), &
                                       netcdf_data_2d, start = (/ 1, 1 /), &
                                       count = (/ mask_size(mid,1), &
                                                  mask_size(mid,2) /) )
               CALL netcdf_handle_error( 'netcdf_define_header', 555 )

               DO  j = 1, mask_size(mid,2)
                  DO  i = 1, mask_size(mid,1)
                     netcdf_data_2d(i,j) = init_model%origin_y                 &
                           - sin_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
                           + cos_ra * ( mask_j_global(mid,j) + shift_y ) * dy
                  ENDDO
               ENDDO
             
               nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
                                       id_var_nutm_mask(k,mid,av), &
                                       netcdf_data_2d, start = (/ 1, 1 /), &
                                       count = (/ mask_size(mid,1), &
                                                  mask_size(mid,2) /) )
               CALL netcdf_handle_error( 'netcdf_define_header', 556 )
             
             ENDDO
             DEALLOCATE( netcdf_data_2d )
          ENDIF
!
!--Â  Â  Â  Â Write lon and lat data.
          ALLOCATE( lat(mask_size(mid,1),mask_size(mid,2)) )
          ALLOCATE( lon(mask_size(mid,1),mask_size(mid,2)) )
          cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
          sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )

          DO  k = 0, 2
!               
!--Â  Â  Â  Â  Â  Scalar grid points
             IF ( k == 0 )  THEN
                shift_x = 0.5 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  u grid points
             ELSEIF ( k == 1 )  THEN
                shift_x = 0.0 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  v grid points
             ELSEIF ( k == 2 )  THEN
                shift_x = 0.5 ; shift_y = 0.0
             ENDIF

             DO  j = 1, mask_size(mid,2)
                DO  i = 1, mask_size(mid,1)
                   eutm = init_model%origin_x                               &
                        + cos_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
                        + sin_ra * ( mask_j_global(mid,j) + shift_y ) * dy
                   nutm = init_model%origin_y                               &
                        - sin_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
                        + cos_ra * ( mask_j_global(mid,j) + shift_y ) * dy

                   CALL  convert_utm_to_geographic( crs_list,          &
                                                    eutm, nutm,        &
                                                    lon(i,j), lat(i,j) )
                ENDDO
             ENDDO

             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),           &
                                     id_var_lon_mask(k,mid,av),     &
                                     lon, start = (/ 1, 1 /),       &
                                     count = (/ mask_size(mid,1),   &
                                                mask_size(mid,2) /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 556 )

             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),           &
                                     id_var_lat_mask(k,mid,av),     &
                                     lat, start = (/ 1, 1 /),       &
                                     count = (/ mask_size(mid,1),   &
                                                mask_size(mid,2) /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 556 )
          ENDDO

          DEALLOCATE( lat )
          DEALLOCATE( lon )
!
!--Â  Â  Â  Â Write zu and zw data (vertical axes)
          ALLOCATE( netcdf_data(mask_size(mid,3)) )

          IF ( mask_surface(mid) )  THEN

             netcdf_data = 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 netcdf_handle_error( 'netcdf_define_header', 503 )

             netcdf_data = 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 netcdf_handle_error( 'netcdf_define_header', 504 )

          ELSE

             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 netcdf_handle_error( 'netcdf_define_header', 503 )

             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 netcdf_handle_error( 'netcdf_define_header', 504 )

          ENDIF

          DEALLOCATE( netcdf_data )

!
!--Â  Â  Â  Â In case of non-flat topography write height information
          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
               netcdf_data_format > 4 )  THEN

             ALLOCATE( netcdf_data_2d(mask_size_l(mid,1),mask_size_l(mid,2)) )
             netcdf_data_2d = zu_s_inner( mask_i(mid,:mask_size_l(mid,1)),     &
                                          mask_j(mid,:mask_size_l(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_l(mid,1),            &
                                                mask_size_l(mid,2) /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 505 )

             netcdf_data_2d = zw_w_inner( mask_i(mid,:mask_size_l(mid,1)),     &
                                          mask_j(mid,:mask_size_l(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_l(mid,1),            &
                                                mask_size_l(mid,2) /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 506 )

             DEALLOCATE( netcdf_data_2d )

          ENDIF

          IF ( land_surface )  THEN
!
!--Â  Â  Â  Â  Â  Write zs data (vertical axes for soil model), use negative values
!--Â  Â  Â  Â  Â  to indicate soil depth
             ALLOCATE( netcdf_data(mask_size(mid,3)) )

             netcdf_data = zs( mask_k_global(mid,:mask_size(mid,3)) )

             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zs_mask(mid,av), &
                                     netcdf_data, start = (/ 1 /), &
                                     count = (/ mask_size(mid,3) /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 538 )

             DEALLOCATE( netcdf_data )

          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 <= 200+max_masks )  THEN
             mid = file_id - 200
             av = 0
          ELSE
             mid = file_id - (200+max_masks)
             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 netcdf_handle_error( 'netcdf_define_header', 507 )

          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', 'PA0335', 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_3d', &
                                    id_var_zu_mask(mid,av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 508 )

          nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av),     &
                                           id_var_zu_mask(mid,av),  &
                                           dimids = id_dim_zu_mask_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 509 )
          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 netcdf_handle_error( 'netcdf_define_header', 510 )

          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', 'PA0336', 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 netcdf_handle_error( 'netcdf_define_header', 511 )

          nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av),                &
                                           id_var_time_mask(mid,av),           &
                                           dimids = id_dim_time_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 512 )
          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 netcdf_handle_error( 'netcdf_define_header', 513 )

          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 netcdf_handle_error( 'netcdf_define_header', 514 )

          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', 'PA0337', 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 netcdf_handle_error( 'netcdf_define_header', 515 )
             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 netcdf_handle_error( 'netcdf_define_header', 516 )
          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'title',   &
                                  TRIM( run_description_header ) //            &
                                  TRIM( time_average_text ) )
          CALL netcdf_handle_error( 'netcdf_define_header', 517 )
          nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 518 )
          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', 'PA0338', 0, 0, 0, 6, 0 )
!
!--Â  Â  Â  Â restore original parameter file_id (=formal parameter av) into av
          av = file_id


       CASE ( '3d_new' )

!
!--Â  Â  Â  Â Define some global attributes of the dataset
          IF ( av == 0 )  THEN
             CALL netcdf_create_global_atts( id_set_3d(av), '3d', TRIM( run_description_header ), 62 )
             time_average_text = ' '
          ELSE
             CALL netcdf_create_global_atts( id_set_3d(av), '3d_av', TRIM( run_description_header ), 62 )
             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
             nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'time_avg',   &
                                     TRIM( time_average_text ) )
             CALL netcdf_handle_error( 'netcdf_define_header', 63 )
          ENDIF

!
!--Â  Â  Â  Â Define time coordinate for volume data.
!--Â  Â  Â  Â For parallel output the time dimensions has to be limited, otherwise
!--Â  Â  Â  Â the performance drops significantly.
          IF ( netcdf_data_format < 5 )  THEN
             CALL netcdf_create_dim( id_set_3d(av), 'time', NF90_UNLIMITED,    &
                                     id_dim_time_3d(av), 64 )
          ELSE
             CALL netcdf_create_dim( id_set_3d(av), 'time', ntdim_3d(av),      &
                                     id_dim_time_3d(av), 523 )
          ENDIF

          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_time_3d(av) /),     &
                                  'time', NF90_DOUBLE, id_var_time_3d(av),     &
                                  'seconds since '//TRIM(init_model%origin_time), 'time', 65, 66, 00 )
          CALL netcdf_create_att( id_set_3d(av), id_var_time_3d(av), 'standard_name', 'time', 000)
          CALL netcdf_create_att( id_set_3d(av), id_var_time_3d(av), 'axis', 'T', 000)
!
!--Â  Â  Â  Â Define spatial dimensions and coordinates:
!--Â  Â  Â  Â Define vertical coordinate grid (zu grid)
          CALL netcdf_create_dim( id_set_3d(av), 'zu_3d', nz_do3d-nzb+1,       &
                                  id_dim_zu_3d(av), 67 )
          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zu_3d(av) /),       &
                                  'zu_3d', NF90_DOUBLE, id_var_zu_3d(av),      &
                                  'meters', '', 68, 69, 00 )
!
!--Â  Â  Â  Â Define vertical coordinate grid (zw grid)
          CALL netcdf_create_dim( id_set_3d(av), 'zw_3d', nz_do3d-nzb+1,       &
                                  id_dim_zw_3d(av), 70 )
          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zw_3d(av) /),       &
                                  'zw_3d', NF90_DOUBLE, id_var_zw_3d(av),      &
                                  'meters', '', 71, 72, 00 )
!
!--Â  Â  Â  Â Define x-axis (for scalar position)
          CALL netcdf_create_dim( id_set_3d(av), 'x', nx+1, id_dim_x_3d(av),   &
                                  73 )
          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av) /), 'x',   &
                                  NF90_DOUBLE, id_var_x_3d(av), 'meters', '',  &
                                  74, 75, 00 )
!
!--Â  Â  Â  Â Define x-axis (for u position)
          CALL netcdf_create_dim( id_set_3d(av), 'xu', nx+1, id_dim_xu_3d(av), &
                                  358 )
          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_xu_3d(av) /), 'xu', &
                                  NF90_DOUBLE, id_var_xu_3d(av), 'meters', '', &
                                  359, 360, 000 )
!
!--Â  Â  Â  Â Define y-axis (for scalar position)
          CALL netcdf_create_dim( id_set_3d(av), 'y', ny+1, id_dim_y_3d(av),   &
                                  76 )
          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_y_3d(av) /), 'y',   &
                                  NF90_DOUBLE, id_var_y_3d(av), 'meters', '',  &
                                  77, 78, 00 )
!
!--Â  Â  Â  Â Define y-axis (for v position)
          CALL netcdf_create_dim( id_set_3d(av), 'yv', ny+1, id_dim_yv_3d(av), &
                                  361 )
          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_yv_3d(av) /), 'yv', &
                                  NF90_DOUBLE, id_var_yv_3d(av), 'meters', '', &
                                  362, 363, 000 )
!
!--Â  Â  Â  Â Define UTM and geographic coordinates
          CALL define_geo_coordinates( id_set_3d(av),         &
                  (/ id_dim_x_3d(av), id_dim_xu_3d(av) /),    &
                  (/ id_dim_y_3d(av), id_dim_yv_3d(av) /),    &
                  id_var_eutm_3d(:,av), id_var_nutm_3d(:,av), &
                  id_var_lat_3d(:,av), id_var_lon_3d(:,av)    )
!
!--Â  Â  Â  Â Define coordinate-reference system
          CALL netcdf_create_crs( id_set_3d(av), 000 )
!
!--Â  Â  Â  Â In case of non-flat topography define 2d-arrays containing the height
!--Â  Â  Â  Â information. Only output 2d topography information in case of parallel
!--Â  Â  Â  Â output.
          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
               netcdf_data_format > 4 )  THEN
!
!--Â  Â  Â  Â  Â  Define zusi = zu(nzb_s_inner)
             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av),        &
                                     id_dim_y_3d(av) /), 'zusi', NF90_DOUBLE,  &
                                     id_var_zusi_3d(av), 'meters',             &
                                     'zu(nzb_s_inner)', 413, 414, 415 )
!             
!--Â  Â  Â  Â  Â  Define zwwi = zw(nzb_w_inner)
             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av),        &
                                     id_dim_y_3d(av) /), 'zwwi', NF90_DOUBLE,  &
                                     id_var_zwwi_3d(av), 'meters',             &
                                     'zw(nzb_w_inner)', 416, 417, 418 )

          ENDIF             

          IF ( land_surface )  THEN
!
!--Â  Â  Â  Â  Â  Define vertical coordinate grid (zs grid)
             CALL netcdf_create_dim( id_set_3d(av), 'zs_3d',                   &
                                     nzt_soil-nzb_soil+1, id_dim_zs_3d(av), 70 )
             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zs_3d(av) /),    &
                                     'zs_3d', NF90_DOUBLE, id_var_zs_3d(av),   &
                                     'meters', '', 71, 72, 00 )

          ENDIF

!
!--Â  Â  Â  Â Define the variables
          var_list = ';'
          i = 1

          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
!
!--Â  Â  Â  Â  Â  Temporary solution to account for data output within the new urban 
!--Â  Â  Â  Â  Â  surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
             trimvar = TRIM( do3d(av,i) )
             IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
                trimvar = 'usm_output'
             ENDIF
!
!--Â  Â  Â  Â  Â  Check for the grid
             found = .FALSE.
             SELECT CASE ( trimvar )
!
!--Â  Â  Â  Â  Â  Â  Â Most variables are defined on the scalar grid
                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',   &
                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
                       's', 'theta', 'thetal', 'thetav' )

                   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'

!             
!--             Block of urban surface model outputs   
                CASE ( 'usm_output' )
                   CALL usm_define_netcdf_grid( do3d(av,i), found, &
                                                   grid_x, grid_y, grid_z )

                CASE DEFAULT

                   CALL tcm_define_netcdf_grid( do3d(av,i), found, &
                                                   grid_x, grid_y, grid_z )

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Check for land surface quantities
                   IF ( .NOT. found .AND. land_surface )  THEN
                      CALL lsm_define_netcdf_grid( do3d(av,i), found, grid_x,  &
                                                   grid_y, grid_z )
                   ENDIF
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Check for ocean quantities
                   IF ( .NOT. found  .AND.  ocean_mode )  THEN
                      CALL ocean_define_netcdf_grid( do3d(av,i), found,  &
                                                     grid_x, grid_y, grid_z )
                   ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Check for plant canopy quantities
                   IF ( .NOT. found  .AND.  plant_canopy )  THEN
                      CALL pcm_define_netcdf_grid( do3d(av,i), found, grid_x,  &
                                                   grid_y, grid_z )
                   ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Check for radiation quantities
                   IF ( .NOT. found  .AND.  radiation )  THEN
                      CALL radiation_define_netcdf_grid( do3d(av,i), found,    &
                                                         grid_x, grid_y,       &
                                                         grid_z )
                   ENDIF

!--Â  Â  Â  Â  Â  Â  Â  Â  Check for gust module quantities
                   IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
                      CALL gust_define_netcdf_grid( do3d(av,i), found, grid_x, &
                                                    grid_y, grid_z )
                   ENDIF
!
!--                Check for indoor model quantities                   
                   IF ( .NOT. found .AND. indoor_model ) THEN 
                      CALL im_define_netcdf_grid( do3d(av,i), found,           &
                                                  grid_x, grid_y, grid_z )
                   ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Check for biometeorology quantities
                   IF ( .NOT. found  .AND.  biometeorology )  THEN
                      CALL bio_define_netcdf_grid( do3d(av,i), found,          &
                                                   grid_x, grid_y, grid_z )
                   ENDIF

!
!--                Check for chemistry quantities                   
                   IF ( .NOT. found  .AND.  air_chemistry )  THEN
                      CALL chem_define_netcdf_grid( do3d(av,i), found,         &
                                                    grid_x, grid_y, grid_z )
                   ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Check for SALSA quantities
                   IF ( .NOT. found  .AND.  salsa )  THEN
                      CALL salsa_define_netcdf_grid( do3d(av,i), found, grid_x,&
                                                     grid_y, grid_z )
                   ENDIF                 
!                   
!--Â  Â  Â  Â  Â  Â  Â  Â  Check for user-defined quantities
                   IF ( .NOT. found  .AND.  user_module_enabled )  THEN
                      CALL user_define_netcdf_grid( do3d(av,i), found, grid_x, &
                                                    grid_y, grid_z )
                   ENDIF
                                                 
                   IF ( .NOT. found )  THEN
                      WRITE ( message_string, * ) 'no grid defined for varia', &
                                                  'ble ', TRIM( do3d(av,i) )
                      CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, &
                                    6, 0 )
                   ENDIF

             END SELECT

!
!--Â  Â  Â  Â  Â  Select the respective dimension ids
             IF ( grid_x == 'x' )  THEN
                id_x = id_dim_x_3d(av)
             ELSEIF ( grid_x == 'xu' )  THEN
                id_x = id_dim_xu_3d(av)
             ENDIF

             IF ( grid_y == 'y' )  THEN
                id_y = id_dim_y_3d(av)
             ELSEIF ( grid_y == 'yv' )  THEN
                id_y = id_dim_yv_3d(av)
             ENDIF

             IF ( grid_z == 'zu' )  THEN
                id_z = id_dim_zu_3d(av)
             ELSEIF ( grid_z == 'zw' )  THEN
                id_z = id_dim_zw_3d(av)
             ELSEIF ( grid_z == 'zs' )  THEN
                id_z = id_dim_zs_3d(av)
             ENDIF

!
!--Â  Â  Â  Â  Â  Define the grid
             CALL netcdf_create_var( id_set_3d(av),(/ id_x, id_y, id_z,        &
                                     id_dim_time_3d(av) /), do3d(av,i),        &
                                     nc_precision(4), id_var_do3d(av,i),       &
                                     TRIM( do3d_unit(av,i) ), do3d(av,i), 79,  &
                                     80, 357, .TRUE. )
#if defined( __netcdf4_parallel )
             IF ( netcdf_data_format > 4 )  THEN
!
!--Â  Â  Â  Â  Â  Â  Â Set no fill for every variable to increase performance.
                nc_stat = NF90_DEF_VAR_FILL( id_set_3d(av),     &
                                             id_var_do3d(av,i), &
                                             1, 0 )
                CALL netcdf_handle_error( 'netcdf_define_header', 532 )
!
!--Â  Â  Â  Â  Â  Â  Â Set collective io operations for parallel io
                nc_stat = NF90_VAR_PAR_ACCESS( id_set_3d(av),     &
                                               id_var_do3d(av,i), &
                                               NF90_COLLECTIVE )
                CALL netcdf_handle_error( 'netcdf_define_header', 445 )
             ENDIF
#endif
             var_list = TRIM( var_list ) // TRIM( do3d(av,i) ) // ';'

             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_3d(av), NF90_GLOBAL, 'VAR_LIST', &
                                  var_list )
          CALL netcdf_handle_error( 'netcdf_define_header', 81 )

!
!--Â  Â  Â  Â Set general no fill, otherwise the performance drops significantly for
!--Â  Â  Â  Â parallel output.
          nc_stat = NF90_SET_FILL( id_set_3d(av), NF90_NOFILL, oldmode )
          CALL netcdf_handle_error( 'netcdf_define_header', 528 )

!
!--Â  Â  Â  Â Leave netCDF define mode
          nc_stat = NF90_ENDDEF( id_set_3d(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 82 )

!
!--Â  Â  Â  Â These data are only written by PE0 for parallel output to increase
!--Â  Â  Â  Â the performance.
          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
!
!--Â  Â  Â  Â  Â  Write data for x (shifted by +dx/2) and xu axis
             ALLOCATE( netcdf_data(0:nx) )

             DO  i = 0, nx
                netcdf_data(i) = ( i + 0.5 ) * dx
             ENDDO

             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_x_3d(av),  &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ nx+1 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 83 )

             DO  i = 0, nx
                netcdf_data(i) = i * dx
             ENDDO

             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_xu_3d(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ nx+1 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 385 )

             DEALLOCATE( netcdf_data )

!
!--Â  Â  Â  Â  Â  Write data for y (shifted by +dy/2) and yv axis
             ALLOCATE( netcdf_data(0:ny) )

             DO  i = 0, ny
                netcdf_data(i) = ( i + 0.5_wp ) * dy
             ENDDO

             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_y_3d(av),  &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ ny+1 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 84 )

             DO  i = 0, ny
                netcdf_data(i) = i * dy
             ENDDO

             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_yv_3d(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ ny+1 /))
             CALL netcdf_handle_error( 'netcdf_define_header', 387 )

             DEALLOCATE( netcdf_data )

!
!--Â  Â  Â  Â  Â  Write UTM coordinates
             IF ( init_model%rotation_angle == 0.0_wp )  THEN
!
!--Â  Â  Â  Â  Â  Â  Â 1D in case of no rotation
                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
!
!--Â  Â  Â  Â  Â  Â  Â x coordinates
                ALLOCATE( netcdf_data(0:nx) )
                DO  k = 0, 2
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  Scalar grid points
                   IF ( k == 0 )  THEN
                      shift_x = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  u grid points
                   ELSEIF ( k == 1 )  THEN
                      shift_x = 0.0
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  v grid points
                   ELSEIF ( k == 2 )  THEN
                      shift_x = 0.5
                   ENDIF
                
                   DO  i = 0, nx
                     netcdf_data(i) = init_model%origin_x            &
                                    + cos_ra * ( i + shift_x ) * dx
                   ENDDO
                
                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_eutm_3d(k,av),&
                                           netcdf_data, start = (/ 1 /),   &
                                           count = (/ nx+1 /) )
                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )

                ENDDO
                DEALLOCATE( netcdf_data )
!
!--Â  Â  Â  Â  Â  Â  Â y coordinates
                ALLOCATE( netcdf_data(0:ny) )
                DO  k = 0, 2
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Scalar grid points
                   IF ( k == 0 )  THEN
                      shift_y = 0.5
!
!--Â  Â  Â  Â  Â  Â  Â  Â  u grid points
                   ELSEIF ( k == 1 )  THEN
                      shift_y = 0.5
!
!--Â  Â  Â  Â  Â  Â  Â  Â  v grid points
                   ELSEIF ( k == 2 )  THEN
                      shift_y = 0.0
                   ENDIF

                   DO  j = 0, ny
                      netcdf_data(j) = init_model%origin_y            &
                                     + cos_ra * ( j + shift_y ) * dy
                   ENDDO

                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_nutm_3d(k,av),&
                                           netcdf_data, start = (/ 1 /),   &
                                           count = (/ ny+1 /) )
                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )

                ENDDO
                DEALLOCATE( netcdf_data )

             ELSE
!
!--Â  Â  Â  Â  Â  Â  Â 2D in case of rotation
                ALLOCATE( netcdf_data_2d(0:nx,0:ny) )
                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
                
                DO  k = 0, 2
!               
!--Â  Â  Â  Â  Â  Â  Â  Â Scalar grid points
                  IF ( k == 0 )  THEN
                     shift_x = 0.5 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â  Â u grid points
                  ELSEIF ( k == 1 )  THEN
                     shift_x = 0.0 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â  Â v grid points
                  ELSEIF ( k == 2 )  THEN
                     shift_x = 0.5 ; shift_y = 0.0
                  ENDIF
                
                  DO  j = 0, ny
                     DO  i = 0, nx
                        netcdf_data_2d(i,j) = init_model%origin_x            &
                                            + cos_ra * ( i + shift_x ) * dx  &
                                            + sin_ra * ( j + shift_y ) * dy
                     ENDDO
                  ENDDO
                
                  nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_eutm_3d(k,av),  &
                                          netcdf_data_2d, start = (/ 1, 1 /),   &
                                          count = (/ nx+1, ny+1 /) )
                  CALL netcdf_handle_error( 'netcdf_define_header', 555 )
                
                  DO  j = 0, ny
                     DO  i = 0, nx
                        netcdf_data_2d(i,j) = init_model%origin_y            &
                                            - sin_ra * ( i + shift_x ) * dx  &
                                            + cos_ra * ( j + shift_y ) * dy
                     ENDDO
                  ENDDO
                
                  nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_nutm_3d(k,av),  &
                                          netcdf_data_2d, start = (/ 1, 1 /),   &
                                          count = (/ nx+1, ny+1 /) )
                  CALL netcdf_handle_error( 'netcdf_define_header', 556 )
                
                ENDDO
                DEALLOCATE( netcdf_data_2d )
             ENDIF
!
!--Â  Â  Â  Â  Â  Write zu and zw data (vertical axes)
             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zu_3d(av),  &
                                     zu(nzb:nz_do3d), start = (/ 1 /), &
                                     count = (/ nz_do3d-nzb+1 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 85 )


             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zw_3d(av),  &
                                     zw(nzb:nz_do3d), start = (/ 1 /), &
                                     count = (/ nz_do3d-nzb+1 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 86 )

             IF ( land_surface )  THEN
!
!--Â  Â  Â  Â  Â  Â  Â Write zs grid
                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zs_3d(av),  &
                                        - zs(nzb_soil:nzt_soil), start = (/ 1 /), &
                                        count = (/ nzt_soil-nzb_soil+1 /) )
                CALL netcdf_handle_error( 'netcdf_define_header', 86 )
             ENDIF

          ENDIF
!
!--Â  Â  Â  Â Write lon and lat data. Only for parallel output.
          IF ( netcdf_data_format > 4 )  THEN

             ALLOCATE( lat(nxl:nxr,nys:nyn) )
             ALLOCATE( lon(nxl:nxr,nys:nyn) )
             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
             sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )

             DO  k = 0, 2
!               
!--Â  Â  Â  Â  Â  Â  Â Scalar grid points
                IF ( k == 0 )  THEN
                   shift_x = 0.5 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â u grid points
                ELSEIF ( k == 1 )  THEN
                   shift_x = 0.0 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â v grid points
                ELSEIF ( k == 2 )  THEN
                   shift_x = 0.5 ; shift_y = 0.0
                ENDIF

                DO  j = nys, nyn
                   DO  i = nxl, nxr
                      eutm = init_model%origin_x            &
                           + cos_ra * ( i + shift_x ) * dx  &
                           + sin_ra * ( j + shift_y ) * dy
                      nutm = init_model%origin_y            &
                           - sin_ra * ( i + shift_x ) * dx  &
                           + cos_ra * ( j + shift_y ) * dy

                      CALL  convert_utm_to_geographic( crs_list,          &
                                                       eutm, nutm,        &
                                                       lon(i,j), lat(i,j) )
                   ENDDO
                ENDDO

                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_lon_3d(k,av), &
                                     lon, start = (/ nxl+1, nys+1 /),       &
                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
                CALL netcdf_handle_error( 'netcdf_define_header', 556 )

                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_lat_3d(k,av), &
                                     lat, start = (/ nxl+1, nys+1 /),       &
                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
             ENDDO

             DEALLOCATE( lat )
             DEALLOCATE( lon )

          ENDIF
!
!--Â  Â  Â  Â In case of non-flat topography write height information. Only for
!--Â  Â  Â  Â parallel netcdf output.
          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
               netcdf_data_format > 4 )  THEN

!             IF ( nxr == nx  .AND.  nyn /= ny )  THEN
!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av),     &
!                                        zu_s_inner(nxl:nxr+1,nys:nyn),         &
!                                        start = (/ nxl+1, nys+1 /),            &
!                                        count = (/ nxr-nxl+2, nyn-nys+1 /) )
!             ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av),     &
!                                        zu_s_inner(nxl:nxr,nys:nyn+1),         &
!                                        start = (/ nxl+1, nys+1 /),            &
!                                        count = (/ nxr-nxl+1, nyn-nys+2 /) )
!             ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av),     &
!                                        zu_s_inner(nxl:nxr+1,nys:nyn+1),       &
!                                        start = (/ nxl+1, nys+1 /),            &
!                                        count = (/ nxr-nxl+2, nyn-nys+2 /) )
!             ELSE
                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av),     &
                                        zu_s_inner(nxl:nxr,nys:nyn),           &
                                        start = (/ nxl+1, nys+1 /),            &
                                        count = (/ nxr-nxl+1, nyn-nys+1 /) )
!             ENDIF
             CALL netcdf_handle_error( 'netcdf_define_header', 419 )

!             IF ( nxr == nx  .AND.  nyn /= ny )  THEN
!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av),     &
!                                        zw_w_inner(nxl:nxr+1,nys:nyn),         &
!                                        start = (/ nxl+1, nys+1 /),            &
!                                        count = (/ nxr-nxl+2, nyn-nys+1 /) )
!             ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av),     &
!                                        zw_w_inner(nxl:nxr,nys:nyn+1),         &
!                                        start = (/ nxl+1, nys+1 /),            &
!                                        count = (/ nxr-nxl+1, nyn-nys+2 /) )
!             ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av),     &
!                                        zw_w_inner(nxl:nxr+1,nys:nyn+1),       &
!                                        start = (/ nxl+1, nys+1 /),            &
!                                        count = (/ nxr-nxl+2, nyn-nys+2 /) )
!             ELSE
                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av),     &
                                        zw_w_inner(nxl:nxr,nys:nyn),           &
                                        start = (/ nxl+1, nys+1 /),            &
                                        count = (/ nxr-nxl+1, nyn-nys+1 /) )
!             ENDIF
             CALL netcdf_handle_error( 'netcdf_define_header', 420 )

          ENDIF

       CASE ( '3d_ext' )

!
!--Â  Â  Â  Â 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_3d(av), NF90_GLOBAL, 'VAR_LIST', &
                                  var_list_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 87 )

          var_list = ';'
          i = 1
          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
             var_list = TRIM(var_list) // TRIM( do3d(av,i) ) // ';'
             i = i + 1
          ENDDO

          IF ( av == 0 )  THEN
             var = '(3d)'
          ELSE
             var = '(3d_av)'
          ENDIF

          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
             message_string = 'netCDF file for volume data ' //             &
                              TRIM( var ) // ' from previous run found,' // &
                              '&but this file cannot be extended due to' // &
                              ' variable mismatch.' //                      &
                              '&New file is created instead.'
             CALL message( 'define_netcdf_header', 'PA0245', 0, 1, 0, 6, 0 )
             extend = .FALSE.
             RETURN
          ENDIF

!
!--Â  Â  Â  Â Get and compare the number of vertical gridpoints
          nc_stat = NF90_INQ_VARID( id_set_3d(av), 'zu_3d', id_var_zu_3d(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 88 )

          nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_zu_3d(av), &
                                           dimids = id_dim_zu_3d_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 89 )
          id_dim_zu_3d(av) = id_dim_zu_3d_old(1)

          nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_zu_3d(av), &
                                            len = nz_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 90 )

          IF ( nz_do3d-nzb+1 /= nz_old )  THEN
              message_string = 'netCDF file for volume data ' //             &
                               TRIM( var ) // ' from previous run found,' // &
                               '&but this file cannot be extended due to' // &
                               ' mismatch in number of' //                   &
                               ' vertical grid points (nz_do3d).' //         &
                               '&New file is created instead.'
             CALL message( 'define_netcdf_header', 'PA0246', 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 pl3d..count+1.
!--Â  Â  Â  Â The current time must be larger than the last output time
!--Â  Â  Â  Â on the file.
          nc_stat = NF90_INQ_VARID( id_set_3d(av), 'time', id_var_time_3d(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 91 )

          nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_time_3d(av), &
                                           dimids = id_dim_time_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 92 )

          id_dim_time_3d(av) = id_dim_time_old(1)

          nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_time_3d(av), &
                                            len = ntime_count )
          CALL netcdf_handle_error( 'netcdf_define_header', 93 )

!
!--Â  Â  Â  Â For non-parallel output use the last output time level of the netcdf
!--Â  Â  Â  Â file because the time dimension is unlimited. In case of parallel
!--Â  Â  Â  Â output the variable ntime_count could get the value of 9*10E36 because
!--Â  Â  Â  Â the time dimension is limited.
          IF ( netcdf_data_format < 5 ) do3d_time_count(av) = ntime_count

          nc_stat = NF90_GET_VAR( id_set_3d(av), id_var_time_3d(av), &
                                  last_time_coordinate,              &
                                  start = (/ do3d_time_count(av) /), &
                                  count = (/ 1 /) )
          CALL netcdf_handle_error( 'netcdf_define_header', 94 )

          IF ( last_time_coordinate(1) >= simulated_time )  THEN
             message_string = 'netCDF file for volume data ' //             &
                              TRIM( var ) // ' from previous run found,' // &
                              '&but this file cannot be extended becaus' // &
                              'e the current output time' //                &
                              '&is less or equal than the last output t' // &
                              'ime on this file.' //                        &
                              '&New file is created instead.'
             CALL message( 'define_netcdf_header', 'PA0247', 0, 1, 0, 6, 0 )
             do3d_time_count(av) = 0
             extend = .FALSE.
             RETURN
          ENDIF

          IF ( netcdf_data_format > 4 )  THEN
!
!--Â  Â  Â  Â  Â  Check if the needed number of output time levels is increased
!--Â  Â  Â  Â  Â  compared to the number of time levels in the existing file.
             IF ( ntdim_3d(av) > ntime_count )  THEN
                message_string = 'netCDF file for volume data ' // &
                                 TRIM( var ) // ' from previous run found,' // &
                                 '&but this file cannot be extended becaus' // &
                                 'e the number of output time levels has b' // &
                                 'een increased compared to the previous s' // &
                                 'imulation.' //                               &
                                 '&New file is created instead.'
                CALL message( 'define_netcdf_header', 'PA0388', 0, 1, 0, 6, 0 )
                do3d_time_count(av) = 0
                extend = .FALSE.
!
!--Â  Â  Â  Â  Â  Â  Â Recalculate the needed time levels for the new file.
                IF ( av == 0 )  THEN
                   ntdim_3d(0) = CEILING(                               &
                           ( end_time - MAX( skip_time_do3d,            &
                                             simulated_time_at_begin )  &
                           ) / dt_do3d )
                   IF ( do3d_at_begin )  ntdim_3d(0) = ntdim_3d(0) + 1
                ELSE
                   ntdim_3d(1) = CEILING(                               &
                           ( end_time - MAX( skip_time_data_output_av,  &
                                             simulated_time_at_begin )  &
                           ) / dt_data_output_av )
                ENDIF
                RETURN
             ENDIF
          ENDIF

!
!--Â  Â  Â  Â Dataset seems to be extendable.
!--Â  Â  Â  Â Now get the variable ids.
          i = 1
          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
             nc_stat = NF90_INQ_VARID( id_set_3d(av), TRIM( do3d(av,i) ), &
                                       id_var_do3d(av,i) )
             CALL netcdf_handle_error( 'netcdf_define_header', 95 )
#if defined( __netcdf4_parallel )
!
!--Â  Â  Â  Â  Â  Set collective io operations for parallel io
             IF ( netcdf_data_format > 4 )  THEN
                nc_stat = NF90_VAR_PAR_ACCESS( id_set_3d(av),     &
                                               id_var_do3d(av,i), &
                                               NF90_COLLECTIVE )
                CALL netcdf_handle_error( 'netcdf_define_header', 453 )
             ENDIF
#endif
             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. 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_3d(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 429 )
          nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'title', &
                                  TRIM( run_description_header ) //    &
                                  TRIM( time_average_text ) )
          CALL netcdf_handle_error( 'netcdf_define_header', 96 )
          nc_stat = NF90_ENDDEF( id_set_3d(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 430 )
          message_string = 'netCDF file for volume data ' //             &
                           TRIM( var ) // ' from previous run found.' // &
                           '&This file will be extended.'
          CALL message( 'define_netcdf_header', 'PA0248', 0, 0, 0, 6, 0 )


       CASE ( 'ag_new' )

!
!--Â  Â  Â  Â Define some global attributes of the dataset
          nc_stat = NF90_PUT_ATT( id_set_agt, NF90_GLOBAL, 'title', &
                                  TRIM( run_description_header ) )
          CALL netcdf_handle_error( 'netcdf_define_header', 330 )
!
!--Â  Â  Â  Â Switch for unlimited time dimension
          IF ( agent_time_unlimited ) THEN
             CALL netcdf_create_dim( id_set_agt, 'time', NF90_UNLIMITED,       &
                                     id_dim_time_agt, 331 )
          ELSE
             CALL netcdf_create_dim( id_set_agt, 'time',                       &
                                     INT( ( MIN( multi_agent_system_end,       &
                                                 end_time ) -                  &
                                            multi_agent_system_start ) /       &
                                            dt_write_agent_data * 1.1 ),       &
                                     id_dim_time_agt, 331 )
          ENDIF

          CALL netcdf_create_var( id_set_agt, (/ id_dim_time_agt /), 'time',   &
                                  NF90_REAL4, id_var_time_agt, 'seconds since '//TRIM(init_model%origin_time), 'time',  &
                                  332, 333, 000 )
          CALL netcdf_create_att( id_set_agt, id_var_time_agt, 'standard_name', 'time', 000)
          CALL netcdf_create_att( id_set_agt, id_var_time_agt, 'axis', 'T', 000)

          CALL netcdf_create_dim( id_set_agt, 'agent_number',                  &
                                  dim_size_agtnum, id_dim_agtnum, 334 )

          CALL netcdf_create_var( id_set_agt, (/ id_dim_agtnum /),             &
                                  'agent_number', NF90_REAL4,                  &
                                  id_var_agtnum, 'agent number', '', 335,      &
                                  336, 000 )
!
!--Â  Â  Â  Â Define variable which contains the real number of agents in use
          CALL netcdf_create_var( id_set_agt, (/ id_dim_time_agt /),           &
                                  'real_num_of_agt', NF90_REAL4,               &
                                  id_var_rnoa_agt, 'agent number', '', 337,    &
                                  338, 000 )
          i = 1
          CALL netcdf_create_var( id_set_agt, (/ id_dim_agtnum,                &
                                  id_dim_time_agt /), agt_var_names(i),        &
                                  NF90_DOUBLE, id_var_agt(i),                  &
                                  TRIM( agt_var_units(i) ),                    &
                                  TRIM( agt_var_names(i) ), 339, 340, 341 )
!
!--Â  Â  Â  Â Define the variables
          DO  i = 2, 6
             CALL netcdf_create_var( id_set_agt, (/ id_dim_agtnum,             &
                                     id_dim_time_agt /), agt_var_names(i),     &
                                     NF90_REAL4, id_var_agt(i),                &
                                     TRIM( agt_var_units(i) ),                 &
                                     TRIM( agt_var_names(i) ), 339, 340, 341 )

          ENDDO
!
!--Â  Â  Â  Â Define vars for biometeorology
          CALL netcdf_create_var( id_set_agt, (/ id_dim_agtnum,                &
                                  id_dim_time_agt /), agt_var_names(9),        &
                                  nc_precision(8), id_var_agt(9),              &
                                  TRIM( agt_var_units(9) ),                    &
                                  TRIM( agt_var_names(9) ), 339, 340, 341 )

!
!--Â  Â  Â  Â Leave netCDF define mode
          nc_stat = NF90_ENDDEF( id_set_agt )
          CALL netcdf_handle_error( 'netcdf_define_header', 342 )


!        CASE ( 'ag_ext' )
! !+?agent extend output for restart runs has to be adapted
! 
! !
! !--Â  Â  Â  Â Get the id of the time coordinate (unlimited coordinate) and its
! !--Â  Â  Â  Â last index on the file. The next time level is prt..count+1.
! !--Â  Â  Â  Â The current time must be larger than the last output time
! !--Â  Â  Â  Â on the file.
!           nc_stat = NF90_INQ_VARID( id_set_agt, 'time', id_var_time_agt )
!           CALL netcdf_handle_error( 'netcdf_define_header', 343 )
! 
!           nc_stat = NF90_INQUIRE_VARIABLE( id_set_agt, id_var_time_agt, &
!                                            dimids = id_dim_time_old )
!           CALL netcdf_handle_error( 'netcdf_define_header', 344 )
!           id_dim_time_agt = id_dim_time_old(1)
! 
!           nc_stat = NF90_INQUIRE_DIMENSION( id_set_agt, id_dim_time_agt, &
!                                             len = agt_time_count )
!           CALL netcdf_handle_error( 'netcdf_define_header', 345 )
! 
!           nc_stat = NF90_GET_VAR( id_set_agt, id_var_time_agt,  &
!                                   last_time_coordinate,         &
!                                   start = (/ agt_time_count /), &
!                                   count = (/ 1 /) )
!           CALL netcdf_handle_error( 'netcdf_define_header', 346 )
! 
!           IF ( last_time_coordinate(1) >= simulated_time )  THEN
!              message_string = 'netCDF file for agents ' //                  &
!                               'from previous run found,' //                 &
!                               '&but this file cannot be extended becaus' // &
!                               'e the current output time' //                &
!                               '&is less or equal than the last output t' // &
!                               'ime on this file.' //                        &
!                               '&New file is created instead.'
!              CALL message( 'define_netcdf_header', 'PA0265', 0, 1, 0, 6, 0 )
!              agt_time_count = 0
!              extend = .FALSE.
!              RETURN
!           ENDIF
! 
! !
! !--Â  Â  Â  Â Dataset seems to be extendable.
! !--Â  Â  Â  Â Now get the variable ids.
!           nc_stat = NF90_INQ_VARID( id_set_agt, 'real_num_of_agt', &
!                                     id_var_rnoa_agt )
!           CALL netcdf_handle_error( 'netcdf_define_header', 347 )
! 
!           DO  i = 1, 17
! 
!              nc_stat = NF90_INQ_VARID( id_set_agt, agt_var_names(i), &
!                                        id_var_prt(i) )
!              CALL netcdf_handle_error( 'netcdf_define_header', 348 )
! 
!           ENDDO
! 
!           message_string = 'netCDF file for particles ' // &
!                            'from previous run found.' //   &
!                            '&This file will be extended.'
!           CALL message( 'define_netcdf_header', 'PA0266', 0, 0, 0, 6, 0 )
          

       CASE ( 'xy_new' )

!
!--Â  Â  Â  Â Define some global attributes of the dataset
          IF ( av == 0 )  THEN
             CALL netcdf_create_global_atts( id_set_xy(av), 'xy', TRIM( run_description_header ), 97 )
             time_average_text = ' '
          ELSE
             CALL netcdf_create_global_atts( id_set_xy(av), 'xy_av', TRIM( run_description_header ), 97 )
             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
             nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'time_avg',   &
                                     TRIM( time_average_text ) )
             CALL netcdf_handle_error( 'netcdf_define_header', 98 )
          ENDIF

!
!--Â  Â  Â  Â Define time coordinate for xy sections.
!--Â  Â  Â  Â For parallel output the time dimensions has to be limited, otherwise
!--Â  Â  Â  Â the performance drops significantly.
          IF ( netcdf_data_format < 5 )  THEN
             CALL netcdf_create_dim( id_set_xy(av), 'time', NF90_UNLIMITED,    &
                                     id_dim_time_xy(av), 99 )
          ELSE
             CALL netcdf_create_dim( id_set_xy(av), 'time', ntdim_2d_xy(av),   &
                                     id_dim_time_xy(av), 524 )
          ENDIF

          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_time_xy(av) /),     &
                                  'time', NF90_DOUBLE, id_var_time_xy(av),     &
                                  'seconds since '//TRIM(init_model%origin_time), 'time', 100, 101, 000 )
          CALL netcdf_create_att( id_set_xy(av), id_var_time_xy(av), 'standard_name', 'time', 000)
          CALL netcdf_create_att( id_set_xy(av), id_var_time_xy(av), 'axis', 'T', 000)
!
!--Â  Â  Â  Â Define the spatial dimensions and coordinates for xy-sections.
!--Â  Â  Â  Â First, determine the number of horizontal sections to be written.
          IF ( section(1,1) == -9999 )  THEN
             RETURN
          ELSE
             ns = 1
             DO WHILE ( section(ns,1) /= -9999  .AND.  ns <= 100 )
                ns = ns + 1
             ENDDO
             ns = ns - 1
          ENDIF

!
!--Â  Â  Â  Â Define vertical coordinate grid (zu grid)
          CALL netcdf_create_dim( id_set_xy(av), 'zu_xy', ns,                  &
                                  id_dim_zu_xy(av), 102 )
          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu_xy(av) /),       &
                                  'zu_xy', NF90_DOUBLE, id_var_zu_xy(av),      &
                                  'meters', '', 103, 104, 000 )
!
!--Â  Â  Â  Â Define vertical coordinate grid (zw grid)
          CALL netcdf_create_dim( id_set_xy(av), 'zw_xy', ns,                  &
                                  id_dim_zw_xy(av), 105 )
          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zw_xy(av) /),       &
                                  'zw_xy', NF90_DOUBLE, id_var_zw_xy(av),      &
                                  'meters', '', 106, 107, 000 )

          IF ( land_surface )  THEN

             ns_do = 1
             DO WHILE ( section(ns_do,1) /= -9999  .AND.  ns_do < nzs )
                ns_do = ns_do + 1
             ENDDO
!
!--Â  Â  Â  Â  Â  Define vertical coordinate grid (zs grid)
             CALL netcdf_create_dim( id_set_xy(av), 'zs_xy', ns_do,            &
                                     id_dim_zs_xy(av), 539 )
             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zs_xy(av) /),    &
                                     'zs_xy', NF90_DOUBLE, id_var_zs_xy(av),   &
                                     'meters', '', 540, 541, 000 )

          ENDIF

!
!--Â  Â  Â  Â Define a pseudo vertical coordinate grid for the surface variables
!--Â  Â  Â  Â u* and t* to store their height level
          CALL netcdf_create_dim( id_set_xy(av), 'zu1_xy', 1,                  &
                                  id_dim_zu1_xy(av), 108 )
          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu1_xy(av) /),      &
                                  'zu1_xy', NF90_DOUBLE, id_var_zu1_xy(av),    &
                                  'meters', '', 109, 110, 000 )
!
!--Â  Â  Â  Â Define a variable to store the layer indices of the horizontal cross
!--Â  Â  Â  Â sections, too
          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu_xy(av) /),       &
                                  'ind_z_xy', NF90_DOUBLE,                     &
                                  id_var_ind_z_xy(av), 'gridpoints', '', 111,  &
                                  112, 000 )
!
!--Â  Â  Â  Â Define x-axis (for scalar position)
          CALL netcdf_create_dim( id_set_xy(av), 'x', nx+1, id_dim_x_xy(av),   &
                                  113 )
          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av) /), 'x',   &
                                  NF90_DOUBLE, id_var_x_xy(av), 'meters', '',  &
                                  114, 115, 000 )
!
!--Â  Â  Â  Â Define x-axis (for u position)
          CALL netcdf_create_dim( id_set_xy(av), 'xu', nx+1,                   &
                                  id_dim_xu_xy(av), 388 )
          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_xu_xy(av) /), 'xu', &
                                  NF90_DOUBLE, id_var_xu_xy(av), 'meters', '', &
                                  389, 390, 000 )
!
!--Â  Â  Â  Â Define y-axis (for scalar position)
          CALL netcdf_create_dim( id_set_xy(av), 'y', ny+1, id_dim_y_xy(av),   &
                                  116 )
          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_y_xy(av) /), 'y',   &
                                  NF90_DOUBLE, id_var_y_xy(av), 'meters', '',  &
                                  117, 118, 000 )
!
!--Â  Â  Â  Â Define y-axis (for scalar position)
          CALL netcdf_create_dim( id_set_xy(av), 'yv', ny+1,                   &
                                  id_dim_yv_xy(av), 364 )
          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_yv_xy(av) /), 'yv', &
                                  NF90_DOUBLE, id_var_yv_xy(av), 'meters', '', &
                                  365, 366, 000 )
!
!--Â  Â  Â  Â Define UTM and geographic coordinates
          CALL define_geo_coordinates( id_set_xy(av),         &
                  (/ id_dim_x_xy(av), id_dim_xu_xy(av) /),    &
                  (/ id_dim_y_xy(av), id_dim_yv_xy(av) /),    &
                  id_var_eutm_xy(:,av), id_var_nutm_xy(:,av), &
                  id_var_lat_xy(:,av), id_var_lon_xy(:,av)    )
!
!--Â  Â  Â  Â Define coordinate-reference system
          CALL netcdf_create_crs( id_set_xy(av), 000 )
!
!--Â  Â  Â  Â In case of non-flat topography define 2d-arrays containing the height
!--Â  Â  Â  Â information. Only for parallel netcdf output.
          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
               netcdf_data_format > 4  )  THEN
!
!--Â  Â  Â  Â  Â  Define zusi = zu(nzb_s_inner)
             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),        &
                                     id_dim_y_xy(av) /), 'zusi', NF90_DOUBLE,  &
                                     id_var_zusi_xy(av), 'meters',             &
                                     'zu(nzb_s_inner)', 421, 422, 423 )
!             
!--Â  Â  Â  Â  Â  Define zwwi = zw(nzb_w_inner)
             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),        &
                                     id_dim_y_xy(av) /), 'zwwi', NF90_DOUBLE,  &
                                     id_var_zwwi_xy(av), 'meters',             &
                                     'zw(nzb_w_inner)', 424, 425, 426 )

          ENDIF

!
!--Â  Â  Â  Â Define the variables
          var_list = ';'
          i = 1

          DO WHILE ( do2d(av,i)(1:1) /= ' ' )

             IF ( INDEX( do2d(av,i), 'xy' ) /= 0 )  THEN
!
!--Â  Â  Â  Â  Â  Â  Â If there is a star in the variable name (u* or t*), it is a
!--Â  Â  Â  Â  Â  Â  Â surface variable. Define it with id_dim_zu1_xy.
                IF ( INDEX( do2d(av,i), '*' ) /= 0 )  THEN

                   CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),  &
                                           id_dim_y_xy(av), id_dim_zu1_xy(av), &
                                           id_dim_time_xy(av) /), do2d(av,i),  &
                                           nc_precision(1), id_var_do2d(av,i), &
                                           TRIM( do2d_unit(av,i) ),            &
                                           do2d(av,i), 119, 120, 354, .TRUE. )

                ELSE

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Check for the grid
                   found = .FALSE.
                   SELECT CASE ( do2d(av,i) )
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Most variables are defined on the zu grid
                      CASE ( 'e_xy', 'nc_xy', 'nr_xy', 'p_xy',                 &
                             'pc_xy', 'pr_xy', 'prr_xy', 'q_xy',               &
                             'qc_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy',           &
                             'ql_vp_xy', 'qr_xy', 'qv_xy',                     &
                             's_xy',                                           &
                             'theta_xy', 'thetal_xy', 'thetav_xy' )

                         grid_x = 'x'
                         grid_y = 'y'
                         grid_z = 'zu'
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â u grid
                      CASE ( 'u_xy' )

                         grid_x = 'xu'
                         grid_y = 'y'
                         grid_z = 'zu'
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â v grid
                      CASE ( 'v_xy' )

                         grid_x = 'x'
                         grid_y = 'yv'
                         grid_z = 'zu'
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â w grid
                      CASE ( 'w_xy' )

                         grid_x = 'x'
                         grid_y = 'y'
                         grid_z = 'zw'


                      CASE DEFAULT
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Check for land surface quantities
                         IF ( land_surface )  THEN
                            CALL lsm_define_netcdf_grid( do2d(av,i), found,    &
                                                   grid_x, grid_y, grid_z )
                         ENDIF

                         IF ( .NOT. found )  THEN
                            CALL tcm_define_netcdf_grid( do2d(av,i), found,    &
                                                         grid_x, grid_y,       &
                                                         grid_z )
                         ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Check for ocean quantities
                         IF ( .NOT. found  .AND.  ocean_mode )  THEN
                            CALL ocean_define_netcdf_grid( do2d(av,i), found,  &
                                                           grid_x, grid_y,     &
                                                           grid_z )
                         ENDIF
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Check for radiation quantities
                         IF ( .NOT. found  .AND.  radiation )  THEN
                            CALL radiation_define_netcdf_grid( do2d(av,i),     &
                                                         found, grid_x, grid_y,&
                                                         grid_z )
                         ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Check for SALSA quantities
                         IF ( .NOT. found  .AND.  salsa )  THEN
                            CALL salsa_define_netcdf_grid( do2d(av,i), found,  &
                                                           grid_x, grid_y,     &
                                                           grid_z )
                         ENDIF                       

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Check for gust module quantities
                         IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
                            CALL gust_define_netcdf_grid( do2d(av,i), found,   &
                                                          grid_x, grid_y,      &
                                                          grid_z )
                         ENDIF
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Check for biometeorology quantities
                         IF ( .NOT. found  .AND.  biometeorology )  THEN
                            CALL bio_define_netcdf_grid( do2d( av, i), found,  &
                                                         grid_x, grid_y,       &
                                                         grid_z )
                         ENDIF
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Check for chemistry quantities
                         IF ( .NOT. found  .AND.  air_chemistry )  THEN
                            CALL chem_define_netcdf_grid( do2d(av,i), found,   &
                                                          grid_x, grid_y,      &
                                                          grid_z )
                         ENDIF


!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Check for user-defined quantities
                         IF ( .NOT. found  .AND.  user_module_enabled )  THEN
                            CALL user_define_netcdf_grid( do2d(av,i), found,   &
                                                          grid_x, grid_y,      &
                                                          grid_z )
                         ENDIF

                         IF ( .NOT. found )  THEN
                            WRITE ( message_string, * ) 'no grid defined for', &
                                                ' variable ', TRIM( do2d(av,i) )
                            CALL message( 'define_netcdf_header', 'PA0244',    &
                                          0, 1, 0, 6, 0 )
                         ENDIF

                   END SELECT

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Select the respective dimension ids
                   IF ( grid_x == 'x' )  THEN
                      id_x = id_dim_x_xy(av)
                   ELSEIF ( grid_x == 'xu' )  THEN
                      id_x = id_dim_xu_xy(av)
                   ENDIF

                   IF ( grid_y == 'y' )  THEN
                      id_y = id_dim_y_xy(av)
                   ELSEIF ( grid_y == 'yv' )  THEN
                      id_y = id_dim_yv_xy(av)
                   ENDIF

                   IF ( grid_z == 'zu' )  THEN
                      id_z = id_dim_zu_xy(av)
                   ELSEIF ( grid_z == 'zw' )  THEN
                      id_z = id_dim_zw_xy(av)
                   ELSEIF ( grid_z == 'zs' )  THEN
                      id_z = id_dim_zs_xy(av)
                   ELSEIF ( grid_z == 'zu1' )  THEN
                      id_z = id_dim_zu1_xy(av)
                   ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Define the grid
                   CALL netcdf_create_var( id_set_xy(av), (/ id_x, id_y, id_z, &
                                           id_dim_time_xy(av) /), do2d(av,i),  &
                                           nc_precision(1), id_var_do2d(av,i), &
                                           TRIM( do2d_unit(av,i) ),            &
                                           do2d(av,i), 119, 120, 354, .TRUE. )

                ENDIF

#if defined( __netcdf4_parallel )
                IF ( netcdf_data_format > 4 )  THEN
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Set no fill for every variable to increase performance.
                   nc_stat = NF90_DEF_VAR_FILL( id_set_xy(av),     &
                                                id_var_do2d(av,i), &
                                                1, 0 )
                   CALL netcdf_handle_error( 'netcdf_define_header', 533 )
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Set collective io operations for parallel io
                   nc_stat = NF90_VAR_PAR_ACCESS( id_set_xy(av),     &
                                                  id_var_do2d(av,i), &
                                                  NF90_COLLECTIVE )
                   CALL netcdf_handle_error( 'netcdf_define_header', 448 )
                ENDIF
#endif
                var_list = TRIM( var_list) // TRIM( do2d(av,i) ) // ';'

             ENDIF

             i = i + 1

          ENDDO

!
!--Â  Â  Â  Â No arrays to output. Close the netcdf file and return.
          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_xy(av), NF90_GLOBAL, 'VAR_LIST', &
                                  var_list )
          CALL netcdf_handle_error( 'netcdf_define_header', 121 )

!
!--Â  Â  Â  Â Set general no fill, otherwise the performance drops significantly for
!--Â  Â  Â  Â parallel output.
          nc_stat = NF90_SET_FILL( id_set_xy(av), NF90_NOFILL, oldmode )
          CALL netcdf_handle_error( 'netcdf_define_header', 529 )

!
!--Â  Â  Â  Â Leave netCDF define mode
          nc_stat = NF90_ENDDEF( id_set_xy(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 122 )

!
!--Â  Â  Â  Â These data are only written by PE0 for parallel output to increase
!--Â  Â  Â  Â the performance.
          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN

!
!--Â  Â  Â  Â  Â  Write axis data: z_xy, x, y
             ALLOCATE( netcdf_data(1:ns) )

!
!--Â  Â  Â  Â  Â  Write zu data
             DO  i = 1, ns
                IF( section(i,1) == -1 )  THEN
                   netcdf_data(i) = -1.0_wp  ! section averaged along z
                ELSE
                   netcdf_data(i) = zu( section(i,1) )
                ENDIF
             ENDDO
             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu_xy(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ ns /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 123 )

!
!--Â  Â  Â  Â  Â  Write zw data
             DO  i = 1, ns
                IF( section(i,1) == -1 )  THEN
                   netcdf_data(i) = -1.0_wp  ! section averaged along z
                ELSE
                   netcdf_data(i) = zw( section(i,1) )
                ENDIF
             ENDDO
             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zw_xy(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ ns /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 124 )

!
!--Â  Â  Â  Â  Â  Write zs data
             IF ( land_surface )  THEN
                ns_do = 0
                DO  i = 1, ns
                   IF( section(i,1) == -1 )  THEN
                      netcdf_data(i) = 1.0_wp  ! section averaged along z
                      ns_do = ns_do + 1
                   ELSEIF ( section(i,1) < nzs )  THEN
                      netcdf_data(i) = - zs( section(i,1) )
                      ns_do = ns_do + 1
                   ENDIF
                ENDDO

                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zs_xy(av), &
                                        netcdf_data(1:ns_do), start = (/ 1 /),    &
                                        count = (/ ns_do /) )
                CALL netcdf_handle_error( 'netcdf_define_header', 124 )

             ENDIF

!
!--Â  Â  Â  Â  Â  Write gridpoint number data
             netcdf_data(1:ns) = section(1:ns,1)
             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_ind_z_xy(av), &
                                     netcdf_data, start = (/ 1 /),       &
                                     count = (/ ns /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 125 )

             DEALLOCATE( netcdf_data )

!
!--Â  Â  Â  Â  Â  Write the cross section height u*, t*
             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu1_xy(av), &
                                     (/ zu(nzb+1) /), start = (/ 1 /), &
                                     count = (/ 1 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 126 )

!
!--Â  Â  Â  Â  Â  Write data for x (shifted by +dx/2) and xu axis
             ALLOCATE( netcdf_data(0:nx) )

             DO  i = 0, nx
                netcdf_data(i) = ( i + 0.5_wp ) * dx
             ENDDO

             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_x_xy(av), &
                                     netcdf_data, start = (/ 1 /),   &
                                     count = (/ nx+1 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 127 )

             DO  i = 0, nx
                netcdf_data(i) = i * dx
             ENDDO

             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_xu_xy(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ nx+1 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 367 )

             DEALLOCATE( netcdf_data )

!
!--Â  Â  Â  Â  Â  Write data for y (shifted by +dy/2) and yv axis
             ALLOCATE( netcdf_data(0:ny+1) )

             DO  i = 0, ny
                netcdf_data(i) = ( i + 0.5_wp ) * dy
             ENDDO

             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_y_xy(av), &
                                     netcdf_data, start = (/ 1 /),   &
                                     count = (/ ny+1 /))
             CALL netcdf_handle_error( 'netcdf_define_header', 128 )

             DO  i = 0, ny
                netcdf_data(i) = i * dy
             ENDDO

             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_yv_xy(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ ny+1 /))
             CALL netcdf_handle_error( 'netcdf_define_header', 368 )

             DEALLOCATE( netcdf_data )
!
!--Â  Â  Â  Â  Â  Write UTM coordinates
             IF ( init_model%rotation_angle == 0.0_wp )  THEN
!
!--Â  Â  Â  Â  Â  Â  Â 1D in case of no rotation
                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
!
!--Â  Â  Â  Â  Â  Â  Â x coordinates
                ALLOCATE( netcdf_data(0:nx) )
                DO  k = 0, 2
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  Scalar grid points
                   IF ( k == 0 )  THEN
                      shift_x = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  u grid points
                   ELSEIF ( k == 1 )  THEN
                      shift_x = 0.0
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  v grid points
                   ELSEIF ( k == 2 )  THEN
                      shift_x = 0.5
                   ENDIF
                
                   DO  i = 0, nx
                     netcdf_data(i) = init_model%origin_x            &
                                    + cos_ra * ( i + shift_x ) * dx
                   ENDDO
                
                   nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_eutm_xy(k,av),&
                                           netcdf_data, start = (/ 1 /),   &
                                           count = (/ nx+1 /) )
                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )

                ENDDO
                DEALLOCATE( netcdf_data )
!
!--Â  Â  Â  Â  Â  Â  Â y coordinates
                ALLOCATE( netcdf_data(0:ny) )
                DO  k = 0, 2
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Scalar grid points
                   IF ( k == 0 )  THEN
                      shift_y = 0.5
!
!--Â  Â  Â  Â  Â  Â  Â  Â  u grid points
                   ELSEIF ( k == 1 )  THEN
                      shift_y = 0.5
!
!--Â  Â  Â  Â  Â  Â  Â  Â  v grid points
                   ELSEIF ( k == 2 )  THEN
                      shift_y = 0.0
                   ENDIF

                   DO  j = 0, ny
                      netcdf_data(j) = init_model%origin_y            &
                                     + cos_ra * ( j + shift_y ) * dy
                   ENDDO

                   nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_nutm_xy(k,av),&
                                           netcdf_data, start = (/ 1 /),   &
                                           count = (/ ny+1 /) )
                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )

                ENDDO
                DEALLOCATE( netcdf_data )

             ELSE
!
!--Â  Â  Â  Â  Â  Â  Â 2D in case of rotation
                ALLOCATE( netcdf_data_2d(0:nx,0:ny) )
                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
                
                DO  k = 0, 2
!               
!--Â  Â  Â  Â  Â  Â  Â  Â Scalar grid points
                  IF ( k == 0 )  THEN
                     shift_x = 0.5 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â  Â u grid points
                  ELSEIF ( k == 1 )  THEN
                     shift_x = 0.0 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â  Â v grid points
                  ELSEIF ( k == 2 )  THEN
                     shift_x = 0.5 ; shift_y = 0.0
                  ENDIF
                
                  DO  j = 0, ny
                     DO  i = 0, nx
                        netcdf_data_2d(i,j) = init_model%origin_x            &
                                            + cos_ra * ( i + shift_x ) * dx  &
                                            + sin_ra * ( j + shift_y ) * dy
                     ENDDO
                  ENDDO
                
                  nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_eutm_xy(k,av),  &
                                          netcdf_data_2d, start = (/ 1, 1 /),   &
                                          count = (/ nx+1, ny+1 /) )
                  CALL netcdf_handle_error( 'netcdf_define_header', 555 )
                
                  DO  j = 0, ny
                     DO  i = 0, nx
                        netcdf_data_2d(i,j) = init_model%origin_y            &
                                            - sin_ra * ( i + shift_x ) * dx  &
                                            + cos_ra * ( j + shift_y ) * dy
                     ENDDO
                  ENDDO
                
                  nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_nutm_xy(k,av),  &
                                          netcdf_data_2d, start = (/ 1, 1 /),   &
                                          count = (/ nx+1, ny+1 /) )
                  CALL netcdf_handle_error( 'netcdf_define_header', 556 )

                ENDDO
                DEALLOCATE( netcdf_data_2d )
             ENDIF

          ENDIF
!
!--Â  Â  Â  Â Write lon and lat data. Only for parallel output.
          IF ( netcdf_data_format > 4 )  THEN

             ALLOCATE( lat(nxl:nxr,nys:nyn) )
             ALLOCATE( lon(nxl:nxr,nys:nyn) )
             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
             sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )

             DO  k = 0, 2
!               
!--Â  Â  Â  Â  Â  Â  Â Scalar grid points
                IF ( k == 0 )  THEN
                   shift_x = 0.5 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â u grid points
                ELSEIF ( k == 1 )  THEN
                   shift_x = 0.0 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â v grid points
                ELSEIF ( k == 2 )  THEN
                   shift_x = 0.5 ; shift_y = 0.0
                ENDIF

                DO  j = nys, nyn
                   DO  i = nxl, nxr
                      eutm = init_model%origin_x            &
                           + cos_ra * ( i + shift_x ) * dx  &
                           + sin_ra * ( j + shift_y ) * dy
                      nutm = init_model%origin_y            &
                           - sin_ra * ( i + shift_x ) * dx  &
                           + cos_ra * ( j + shift_y ) * dy

                      CALL  convert_utm_to_geographic( crs_list,          &
                                                       eutm, nutm,        &
                                                       lon(i,j), lat(i,j) )
                   ENDDO
                ENDDO

                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_lon_xy(k,av), &
                                     lon, start = (/ nxl+1, nys+1 /),       &
                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
                CALL netcdf_handle_error( 'netcdf_define_header', 556 )

                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_lat_xy(k,av), &
                                     lat, start = (/ nxl+1, nys+1 /),       &
                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
             ENDDO

             DEALLOCATE( lat )
             DEALLOCATE( lon )

          ENDIF
!
!--Â  Â  Â  Â In case of non-flat topography write height information. Only for 
!--Â  Â  Â  Â parallel netcdf output.
          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
               netcdf_data_format > 4  )  THEN

!             IF ( nxr == nx  .AND.  nyn /= ny )  THEN
!                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av),     &
!                                        zu_s_inner(nxl:nxr+1,nys:nyn),         &
!                                        start = (/ nxl+1, nys+1 /),            &
!                                        count = (/ nxr-nxl+2, nyn-nys+1 /) )
!             ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
!                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av),     &
!                                        zu_s_inner(nxl:nxr,nys:nyn+1),         &
!                                        start = (/ nxl+1, nys+1 /),            &
!                                        count = (/ nxr-nxl+1, nyn-nys+2 /) )
!             ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
!                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av),     &
!                                        zu_s_inner(nxl:nxr+1,nys:nyn+1),       &
!                                        start = (/ nxl+1, nys+1 /),            &
!                                        count = (/ nxr-nxl+2, nyn-nys+2 /) )
!             ELSE
                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av),     &
                                        zu_s_inner(nxl:nxr,nys:nyn),           &
                                        start = (/ nxl+1, nys+1 /),            &
                                        count = (/ nxr-nxl+1, nyn-nys+1 /) )
!             ENDIF
             CALL netcdf_handle_error( 'netcdf_define_header', 427 )

!             IF ( nxr == nx  .AND.  nyn /= ny )  THEN
!                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av),     &
!                                        zw_w_inner(nxl:nxr+1,nys:nyn),         &
!                                        start = (/ nxl+1, nys+1 /),            &
!                                        count = (/ nxr-nxl+2, nyn-nys+1 /) )
!             ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
!                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av),     &
!                                        zw_w_inner(nxl:nxr,nys:nyn+1),         &
!                                        start = (/ nxl+1, nys+1 /),            &
!                                        count = (/ nxr-nxl+1, nyn-nys+2 /) )
!             ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
!                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av),     &
!                                        zw_w_inner(nxl:nxr+1,nys:nyn+1),       &
!                                        start = (/ nxl+1, nys+1 /),            &
!                                        count = (/ nxr-nxl+2, nyn-nys+2 /) )
!             ELSE
                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av),     &
                                        zw_w_inner(nxl:nxr,nys:nyn),           &
                                        start = (/ nxl+1, nys+1 /),            &
                                        count = (/ nxr-nxl+1, nyn-nys+1 /) )
!             ENDIF
             CALL netcdf_handle_error( 'netcdf_define_header', 428 )

          ENDIF

       CASE ( 'xy_ext' )

!
!--Â  Â  Â  Â 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_xy(av), NF90_GLOBAL, 'VAR_LIST', &
                                  var_list_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 129 )

          var_list = ';'
          i = 1
          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
             IF ( INDEX( do2d(av,i), 'xy' ) /= 0 )  THEN
                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
             ENDIF
             i = i + 1
          ENDDO

          IF ( av == 0 )  THEN
             var = '(xy)'
          ELSE
             var = '(xy_av)'
          ENDIF

          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
             message_string = 'netCDF file for cross-sections ' //           &
                              TRIM( var ) // ' from previous run found,' //  &
                              '&but this file cannot be extended due to' //  &
                              ' variable mismatch.' //                       &
                              '&New file is created instead.'
             CALL message( 'define_netcdf_header', 'PA0249', 0, 1, 0, 6, 0 )
             extend = .FALSE.
             RETURN
          ENDIF

!
!--Â  Â  Â  Â Calculate the number of current sections
          ns = 1
          DO WHILE ( section(ns,1) /= -9999  .AND.  ns <= 100 )
             ns = ns + 1
          ENDDO
          ns = ns - 1

!
!--Â  Â  Â  Â Get and compare the number of horizontal cross sections
          nc_stat = NF90_INQ_VARID( id_set_xy(av), 'zu_xy', id_var_zu_xy(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 130 )

          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_zu_xy(av), &
                                           dimids = id_dim_zu_xy_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 131 )
          id_dim_zu_xy(av) = id_dim_zu_xy_old(1)

          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_zu_xy(av), &
                                            len = ns_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 132 )

          IF ( ns /= ns_old )  THEN
             message_string = 'netCDF file for cross-sections ' //          &
                              TRIM( var ) // ' from previous run found,' // &
                              '&but this file cannot be extended due to' // &
                              ' mismatch in number of' //                   &
                              ' cross sections.' //                         &
                              '&New file is created instead.'
             CALL message( 'define_netcdf_header', 'PA0250', 0, 1, 0, 6, 0 )
             extend = .FALSE.
             RETURN
          ENDIF

!
!--Â  Â  Â  Â Get and compare the heights of the cross sections
          ALLOCATE( netcdf_data(1:ns_old) )

          nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_zu_xy(av), netcdf_data )
          CALL netcdf_handle_error( 'netcdf_define_header', 133 )

          DO  i = 1, ns
             IF ( section(i,1) /= -1 )  THEN
                IF ( zu(section(i,1)) /= netcdf_data(i) )  THEN
                   message_string = 'netCDF file for cross-sections ' //       &
                               TRIM( var ) // ' from previous run found,' //   &
                               ' but this file cannot be extended' //          &
                               ' due to mismatch in cross' //                  &
                               ' section levels.' //                           &
                               ' New file is created instead.'
                   CALL message( 'define_netcdf_header', 'PA0251',             &
                                                                 0, 1, 0, 6, 0 )
                   extend = .FALSE.
                   RETURN
                ENDIF
             ELSE
                IF ( -1.0_wp /= netcdf_data(i) )  THEN
                   message_string = 'netCDF file for cross-sections ' //       &
                               TRIM( var ) // ' from previous run found,' //   &
                               ' but this file cannot be extended' //          &
                               ' due to mismatch in cross' //                  &
                               ' section levels.' //                           &
                               ' New file is created instead.'
                   CALL message( 'define_netcdf_header', 'PA0251',             &
                                                                 0, 1, 0, 6, 0 )
                   extend = .FALSE.
                   RETURN
                ENDIF
             ENDIF
          ENDDO

          DEALLOCATE( netcdf_data )

!
!--Â  Â  Â  Â Get the id of the time coordinate (unlimited coordinate) and its
!--Â  Â  Â  Â last index on the file. The next time level is do2d..count+1.
!--Â  Â  Â  Â The current time must be larger than the last output time
!--Â  Â  Â  Â on the file.
          nc_stat = NF90_INQ_VARID( id_set_xy(av), 'time', id_var_time_xy(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 134 )

          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_time_xy(av), &
                                           dimids = id_dim_time_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 135 )
          id_dim_time_xy(av) = id_dim_time_old(1)

          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_time_xy(av), &
                                            len = ntime_count )
          CALL netcdf_handle_error( 'netcdf_define_header', 136 )

!
!--Â  Â  Â  Â For non-parallel output use the last output time level of the netcdf
!--Â  Â  Â  Â file because the time dimension is unlimited. In case of parallel
!--Â  Â  Â  Â output the variable ntime_count could get the value of 9*10E36 because
!--Â  Â  Â  Â the time dimension is limited.
          IF ( netcdf_data_format < 5 ) do2d_xy_time_count(av) = ntime_count

          nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_time_xy(av),           &
                                  last_time_coordinate,                        &
                                  start = (/ do2d_xy_time_count(av) /),        &
                                  count = (/ 1 /) )
          CALL netcdf_handle_error( 'netcdf_define_header', 137 )

          IF ( last_time_coordinate(1) >= simulated_time )  THEN
             message_string = 'netCDF file for cross sections ' //             &
                              TRIM( var ) // ' from previous run found,' //    &
                              '&but this file cannot be extended becaus' //    &
                              'e the current output time' //                   &
                              '&is less or equal than the last output t' //    &
                              'ime on this file.' //                           &
                              '&New file is created instead.'
             CALL message( 'define_netcdf_header', 'PA0252', 0, 1, 0, 6, 0 )
             do2d_xy_time_count(av) = 0
             extend = .FALSE.
             RETURN
          ENDIF

          IF ( netcdf_data_format > 4 )  THEN
!
!--Â  Â  Â  Â  Â  Check if the needed number of output time levels is increased
!--Â  Â  Â  Â  Â  compared to the number of time levels in the existing file.
             IF ( ntdim_2d_xy(av) > ntime_count )  THEN
                message_string = 'netCDF file for cross sections ' //          &
                                 TRIM( var ) // ' from previous run found,' // &
                                 '&but this file cannot be extended becaus' // &
                                 'e the number of output time levels has b' // &
                                 'een increased compared to the previous s' // &
                                 'imulation.' //                               &
                                 '&New file is created instead.'
                CALL message( 'define_netcdf_header', 'PA0389', 0, 1, 0, 6, 0 )
                do2d_xy_time_count(av) = 0
                extend = .FALSE.
!
!--Â  Â  Â  Â  Â  Â  Â Recalculate the needed time levels for the new file.
                IF ( av == 0 )  THEN
                   ntdim_2d_xy(0) = CEILING(                            &
                           ( end_time - MAX( skip_time_do2d_xy,         &
                                             simulated_time_at_begin )  &
                           ) / dt_do2d_xy )
                   IF ( do2d_at_begin )  ntdim_2d_xy(0) = ntdim_2d_xy(0) + 1
                ELSE
                   ntdim_2d_xy(1) = CEILING(                            &
                           ( end_time - MAX( skip_time_data_output_av,  &
                                             simulated_time_at_begin )  &
                           ) / dt_data_output_av )
                ENDIF
                RETURN
             ENDIF
          ENDIF

!
!--Â  Â  Â  Â Dataset seems to be extendable.
!--Â  Â  Â  Â Now get the variable ids.
          i = 1
          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
             IF ( INDEX( do2d(av,i), 'xy' ) /= 0 )  THEN
                nc_stat = NF90_INQ_VARID( id_set_xy(av), do2d(av,i), &
                                          id_var_do2d(av,i) )
                CALL netcdf_handle_error( 'netcdf_define_header', 138 )
#if defined( __netcdf4_parallel )
!
!--Â  Â  Â  Â  Â  Â  Â Set collective io operations for parallel io
                IF ( netcdf_data_format > 4 )  THEN
                   nc_stat = NF90_VAR_PAR_ACCESS( id_set_xy(av),     &
                                                  id_var_do2d(av,i), &
                                                  NF90_COLLECTIVE )
                   CALL netcdf_handle_error( 'netcdf_define_header', 454 )
                ENDIF
#endif
             ENDIF
             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_xy(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 431 )
          nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'title',         &
                                  TRIM( run_description_header ) //            &
                                  TRIM( time_average_text ) )
          CALL netcdf_handle_error( 'netcdf_define_header', 139 )
          nc_stat = NF90_ENDDEF( id_set_xy(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 432 )
          message_string = 'netCDF file for cross-sections ' //                &
                            TRIM( var ) // ' from previous run found.' //      &
                           '&This file will be extended.'
          CALL message( 'define_netcdf_header', 'PA0253', 0, 0, 0, 6, 0 )
          

       CASE ( 'xz_new' )

!
!--Â  Â  Â  Â Define some global attributes of the dataset
          IF ( av == 0 )  THEN
             CALL netcdf_create_global_atts( id_set_xz(av), 'xz', TRIM( run_description_header ), 140 )
             time_average_text = ' '
          ELSE
             CALL netcdf_create_global_atts( id_set_xz(av), 'xz_av', TRIM( run_description_header ), 140 )
             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
             nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'time_avg',   &
                                     TRIM( time_average_text ) )
             CALL netcdf_handle_error( 'netcdf_define_header', 141 )
          ENDIF

!
!--Â  Â  Â  Â Define time coordinate for xz sections.
!--Â  Â  Â  Â For parallel output the time dimensions has to be limited, otherwise
!--Â  Â  Â  Â the performance drops significantly.
          IF ( netcdf_data_format < 5 )  THEN
             CALL netcdf_create_dim( id_set_xz(av), 'time', NF90_UNLIMITED,    &
                                     id_dim_time_xz(av), 142 )
          ELSE
             CALL netcdf_create_dim( id_set_xz(av), 'time', ntdim_2d_xz(av),   &
                                     id_dim_time_xz(av), 525 )
          ENDIF

          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_time_xz(av) /),     &
                                  'time', NF90_DOUBLE, id_var_time_xz(av),     &
                                  'seconds since '//TRIM(init_model%origin_time), 'time', 143, 144, 000 )
          CALL netcdf_create_att( id_set_xz(av), id_var_time_xz(av), 'standard_name', 'time', 000)
          CALL netcdf_create_att( id_set_xz(av), id_var_time_xz(av), 'axis', 'T', 000)
!
!--Â  Â  Â  Â Define the spatial dimensions and coordinates for xz-sections.
!--Â  Â  Â  Â First, determine the number of vertical sections to be written.
          IF ( section(1,2) == -9999 )  THEN
             RETURN
          ELSE
             ns = 1
             DO WHILE ( section(ns,2) /= -9999  .AND.  ns <= 100 )
                ns = ns + 1
             ENDDO
             ns = ns - 1
          ENDIF

!
!--Â  Â  Â  Â Define y-axis (for scalar position)
          CALL netcdf_create_dim( id_set_xz(av), 'y_xz', ns, id_dim_y_xz(av),  &
                                  145 )
          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_y_xz(av) /),        &
                                  'y_xz', NF90_DOUBLE, id_var_y_xz(av),        &
                                  'meters', '', 146, 147, 000 )
!
!--Â  Â  Â  Â Define y-axis (for v position)
          CALL netcdf_create_dim( id_set_xz(av), 'yv_xz', ns,                  &
                                  id_dim_yv_xz(av), 369 )
          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_yv_xz(av) /),       &
                                  'yv_xz', NF90_DOUBLE, id_var_yv_xz(av),      &
                                  'meters', '', 370, 371, 000 )
!
!--Â  Â  Â  Â Define a variable to store the layer indices of the vertical cross
!--Â  Â  Â  Â sections
          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_y_xz(av) /),        &
                                  'ind_y_xz', NF90_DOUBLE,                     &
                                  id_var_ind_y_xz(av), 'gridpoints', '', 148,  &
                                  149, 000 )
!
!--Â  Â  Â  Â Define x-axis (for scalar position)
          CALL netcdf_create_dim( id_set_xz(av), 'x', nx+1, id_dim_x_xz(av),   &
                                  150 )
          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_x_xz(av) /), 'x',   &
                                  NF90_DOUBLE, id_var_x_xz(av), 'meters', '',  &
                                  151, 152, 000 )
!
!--Â  Â  Â  Â Define x-axis (for u position)
          CALL netcdf_create_dim( id_set_xz(av), 'xu', nx+1, id_dim_xu_xz(av), &
                                  372 )
          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_xu_xz(av) /), 'xu', &
                                  NF90_DOUBLE, id_var_xu_xz(av), 'meters', '', &
                                  373, 374, 000 )
!
!--Â  Â  Â  Â Define the three z-axes (zu, zw, and zs)
          CALL netcdf_create_dim( id_set_xz(av), 'zu', nz+2, id_dim_zu_xz(av), &
                                  153 )
          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zu_xz(av) /), 'zu', &
                                  NF90_DOUBLE, id_var_zu_xz(av), 'meters', '', &
                                  154, 155, 000 )
          CALL netcdf_create_dim( id_set_xz(av), 'zw', nz+2, id_dim_zw_xz(av), &
                                  156 )
          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zw_xz(av) /), 'zw', &
                                  NF90_DOUBLE, id_var_zw_xz(av), 'meters', '', &
                                  157, 158, 000 )
!
!--Â  Â  Â  Â Define UTM and geographic coordinates
          CALL define_geo_coordinates( id_set_xz(av),         &
                  (/ id_dim_x_xz(av), id_dim_xu_xz(av) /),    &
                  (/ id_dim_y_xz(av), id_dim_yv_xz(av) /),    &
                  id_var_eutm_xz(:,av), id_var_nutm_xz(:,av), &
                  id_var_lat_xz(:,av), id_var_lon_xz(:,av)    )
!
!--Â  Â  Â  Â Define coordinate-reference system
          CALL netcdf_create_crs( id_set_xz(av), 000 )

          IF ( land_surface )  THEN

             CALL netcdf_create_dim( id_set_xz(av), 'zs', nzs,                 &
                                     id_dim_zs_xz(av), 542 )
             CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zs_xz(av) /),    &
                                     'zs', NF90_DOUBLE, id_var_zs_xz(av),      &
                                     'meters', '', 543, 544, 000 )

          ENDIF

!
!--Â  Â  Â  Â Define the variables
          var_list = ';'
          i = 1

          DO WHILE ( do2d(av,i)(1:1) /= ' ' )

             IF ( INDEX( do2d(av,i), 'xz' ) /= 0 )  THEN

!
!--Â  Â  Â  Â  Â  Â  Â Check for the grid
                found = .FALSE.
                SELECT CASE ( do2d(av,i) )
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Most variables are defined on the zu grid
                   CASE ( 'e_xz', 'nc_xz', 'nr_xz', 'p_xz', 'pc_xz',           &
                          'pr_xz', 'prr_xz', 'q_xz', 'qc_xz',                  &
                          'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz',  &
                          'qv_xz', 's_xz',                                     &
                          'theta_xz', 'thetal_xz', 'thetav_xz'                 )

                      grid_x = 'x'
                      grid_y = 'y'
                      grid_z = 'zu'
!
!--Â  Â  Â  Â  Â  Â  Â  Â  u grid
                   CASE ( 'u_xz' )

                      grid_x = 'xu'
                      grid_y = 'y'
                      grid_z = 'zu'
!
!--Â  Â  Â  Â  Â  Â  Â  Â  v grid
                   CASE ( 'v_xz' )

                      grid_x = 'x'
                      grid_y = 'yv'
                      grid_z = 'zu'
!
!--Â  Â  Â  Â  Â  Â  Â  Â  w grid
                   CASE ( 'w_xz' )

                      grid_x = 'x'
                      grid_y = 'y'
                      grid_z = 'zw'

                   CASE DEFAULT

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for land surface quantities
                      IF ( land_surface )  THEN
                         CALL lsm_define_netcdf_grid( do2d(av,i), found,       &
                                                      grid_x, grid_y, grid_z )
                      ENDIF

                      IF ( .NOT. found )  THEN
                         CALL tcm_define_netcdf_grid( do2d(av,i), found,       &
                                                      grid_x, grid_y, grid_z )
                      ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for ocean quantities
                      IF ( .NOT. found  .AND.  ocean_mode )  THEN
                         CALL ocean_define_netcdf_grid( do2d(av,i), found,  &
                                                        grid_x, grid_y, grid_z )
                      ENDIF
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for radiation quantities
                      IF ( .NOT. found  .AND.  radiation )  THEN
                         CALL radiation_define_netcdf_grid( do2d(av,i), found, &
                                                            grid_x, grid_y,    &
                                                            grid_z )
                      ENDIF
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for SALSA quantities
                      IF ( .NOT. found  .AND.  salsa )  THEN
                         CALL salsa_define_netcdf_grid( do2d(av,i), found,     &
                                                        grid_x, grid_y, grid_z )
                      ENDIF                         

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for gust module quantities
                      IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
                         CALL gust_define_netcdf_grid( do2d(av,i), found,      &
                                                       grid_x, grid_y, grid_z )
                      ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for chemistry quantities
                      IF ( .NOT. found  .AND.  air_chemistry )  THEN
                         CALL chem_define_netcdf_grid( do2d(av,i), found,      &
                                                       grid_x, grid_y,         &
                                                       grid_z )
                      ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for user-defined quantities
                      IF ( .NOT. found  .AND.  user_module_enabled )  THEN
                         CALL user_define_netcdf_grid( do2d(av,i), found,      &
                                                       grid_x, grid_y, grid_z )
                      ENDIF

                      IF ( .NOT. found )  THEN
                         WRITE ( message_string, * ) 'no grid defined for',    &
                                                ' variable ', TRIM( do2d(av,i) )
                         CALL message( 'define_netcdf_header', 'PA0244',       &
                                       0, 1, 0, 6, 0 )
                      ENDIF

                END SELECT

!
!--Â  Â  Â  Â  Â  Â  Â Select the respective dimension ids
                IF ( grid_x == 'x' )  THEN
                   id_x = id_dim_x_xz(av)
                ELSEIF ( grid_x == 'xu' )  THEN
                   id_x = id_dim_xu_xz(av)
                ENDIF

                IF ( grid_y == 'y' )  THEN
                   id_y = id_dim_y_xz(av)
                ELSEIF ( grid_y == 'yv' )  THEN
                   id_y = id_dim_yv_xz(av)
                ENDIF

                IF ( grid_z == 'zu' )  THEN
                   id_z = id_dim_zu_xz(av)
                ELSEIF ( grid_z == 'zw' )  THEN
                   id_z = id_dim_zw_xz(av)
                ELSEIF ( grid_z == 'zs' )  THEN
                   id_z = id_dim_zs_xz(av)
                ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â Define the grid
                CALL netcdf_create_var( id_set_xz(av), (/ id_x, id_y, id_z,    &
                                        id_dim_time_xz(av) /), do2d(av,i),     &
                                        nc_precision(2), id_var_do2d(av,i),    &
                                        TRIM( do2d_unit(av,i) ), do2d(av,i),   &
                                        159, 160, 355, .TRUE. )

#if defined( __netcdf4_parallel )

                IF ( netcdf_data_format > 4 )  THEN
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Set no fill for every variable to increase performance.
                   nc_stat = NF90_DEF_VAR_FILL( id_set_xz(av),     &
                                                id_var_do2d(av,i), &
                                                1, 0 )
                   CALL netcdf_handle_error( 'netcdf_define_header', 534 )
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Set independent io operations for parallel io. Collective io
!--Â  Â  Â  Â  Â  Â  Â  Â  is only allowed in case of a 1d-decomposition along x,
!--Â  Â  Â  Â  Â  Â  Â  Â  because otherwise, not all PEs have output data.
                   IF ( npey == 1 )  THEN
                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
                                                     id_var_do2d(av,i), &
                                                     NF90_COLLECTIVE )
                   ELSE
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Test simulations showed that the output of cross sections
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â by all PEs in data_output_2d using NF90_COLLECTIVE is 
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â faster than the output by the first row of PEs in
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â x-direction using NF90_INDEPENDENT.
                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),    & 
                                                    id_var_do2d(av,i), &
                                                    NF90_COLLECTIVE )
!                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
!                                                     id_var_do2d(av,i), &
!                                                     NF90_INDEPENDENT )
                   ENDIF
                   CALL netcdf_handle_error( 'netcdf_define_header', 449 )
                ENDIF
#endif
                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'

             ENDIF

             i = i + 1

          ENDDO

!
!--Â  Â  Â  Â No arrays to output. Close the netcdf file and return.
          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_xz(av), NF90_GLOBAL, 'VAR_LIST', &
                                  var_list )
          CALL netcdf_handle_error( 'netcdf_define_header', 161 )

!
!--Â  Â  Â  Â Set general no fill, otherwise the performance drops significantly for
!--Â  Â  Â  Â parallel output.
          nc_stat = NF90_SET_FILL( id_set_xz(av), NF90_NOFILL, oldmode )
          CALL netcdf_handle_error( 'netcdf_define_header', 530 )

!
!--Â  Â  Â  Â Leave netCDF define mode
          nc_stat = NF90_ENDDEF( id_set_xz(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 162 )

!
!--Â  Â  Â  Â These data are only written by PE0 for parallel output to increase
!--Â  Â  Â  Â the performance.
          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN

!
!--Â  Â  Â  Â  Â  Write axis data: y_xz, x, zu, zw
             ALLOCATE( netcdf_data(1:ns) )

!
!--Â  Â  Â  Â  Â  Write y_xz data (shifted by +dy/2)
             DO  i = 1, ns
                IF( section(i,2) == -1 )  THEN
                   netcdf_data(i) = -1.0_wp  ! section averaged along y
                ELSE
                   netcdf_data(i) = ( section(i,2) + 0.5_wp ) * dy
                ENDIF
             ENDDO
             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_y_xz(av), &
                                     netcdf_data, start = (/ 1 /),   &
                                     count = (/ ns /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 163 )

!
!--Â  Â  Â  Â  Â  Write yv_xz data
             DO  i = 1, ns
                IF( section(i,2) == -1 )  THEN
                   netcdf_data(i) = -1.0_wp  ! section averaged along y
                ELSE
                   netcdf_data(i) = section(i,2) * dy
                ENDIF
             ENDDO
             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_yv_xz(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ ns /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 375 )

!
!--Â  Â  Â  Â  Â  Write gridpoint number data
             netcdf_data(1:ns) = section(1:ns,2)
             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_ind_y_xz(av), &
                                     netcdf_data, start = (/ 1 /),       &
                                     count = (/ ns /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 164 )


             DEALLOCATE( netcdf_data )

!
!--Â  Â  Â  Â  Â  Write data for x (shifted by +dx/2) and xu axis
             ALLOCATE( netcdf_data(0:nx) )

             DO  i = 0, nx
                netcdf_data(i) = ( i + 0.5_wp ) * dx
             ENDDO

             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_x_xz(av), &
                                     netcdf_data, start = (/ 1 /),   &
                                     count = (/ nx+1 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 165 )

             DO  i = 0, nx
                netcdf_data(i) = i * dx
             ENDDO

             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_xu_xz(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ nx+1 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 377 )

             DEALLOCATE( netcdf_data )

!
!--Â  Â  Â  Â  Â  Write zu and zw data (vertical axes)
             ALLOCATE( netcdf_data(0:nz+1) )

             netcdf_data(0:nz+1) = zu(nzb:nzt+1)
             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zu_xz(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ nz+2 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 166 )

             netcdf_data(0:nz+1) = zw(nzb:nzt+1)
             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zw_xz(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ nz+2 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 167 )

!
!--Â  Â  Â  Â  Â  Write zs data
             IF ( land_surface )  THEN
                netcdf_data(0:nzs-1) = - zs(nzb_soil:nzt_soil)
                nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zs_xz(av), &
                                        netcdf_data(0:nzs), start = (/ 1 /),    &
                                        count = (/ nzt_soil-nzb_soil+1 /) )
               CALL netcdf_handle_error( 'netcdf_define_header', 548 )
             ENDIF

             DEALLOCATE( netcdf_data )
!
!--Â  Â  Â  Â  Â  Write UTM coordinates
             IF ( init_model%rotation_angle == 0.0_wp )  THEN
!
!--Â  Â  Â  Â  Â  Â  Â 1D in case of no rotation
                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
!
!--Â  Â  Â  Â  Â  Â  Â x coordinates
                ALLOCATE( netcdf_data(0:nx) )
                DO  k = 0, 2
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  Scalar grid points
                   IF ( k == 0 )  THEN
                      shift_x = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  u grid points
                   ELSEIF ( k == 1 )  THEN
                      shift_x = 0.0
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  v grid points
                   ELSEIF ( k == 2 )  THEN
                      shift_x = 0.5
                   ENDIF
                
                   DO  i = 0, nx
                     netcdf_data(i) = init_model%origin_x            &
                                    + cos_ra * ( i + shift_x ) * dx
                   ENDDO
                
                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_eutm_xz(k,av),&
                                           netcdf_data, start = (/ 1 /),   &
                                           count = (/ nx+1 /) )
                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )

                ENDDO
                DEALLOCATE( netcdf_data )
!
!--Â  Â  Â  Â  Â  Â  Â y coordinates
                ALLOCATE( netcdf_data(1:ns) )
                DO  k = 0, 2
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Scalar grid points
                   IF ( k == 0 )  THEN
                      shift_y = 0.5
!
!--Â  Â  Â  Â  Â  Â  Â  Â  u grid points
                   ELSEIF ( k == 1 )  THEN
                      shift_y = 0.5
!
!--Â  Â  Â  Â  Â  Â  Â  Â  v grid points
                   ELSEIF ( k == 2 )  THEN
                      shift_y = 0.0
                   ENDIF

                   DO  i = 1, ns
                      IF( section(i,2) == -1 )  THEN
                         netcdf_data(i) = -1.0_wp  ! section averaged along y
                      ELSE
                         netcdf_data(i) = init_model%origin_y &
                                     + cos_ra * ( section(i,2) + shift_y ) * dy
                      ENDIF
                   ENDDO

                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_nutm_xz(k,av),&
                                           netcdf_data, start = (/ 1 /),   &
                                           count = (/ ns /) )
                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )

                ENDDO
                DEALLOCATE( netcdf_data )

             ELSE
!
!--Â  Â  Â  Â  Â  Â  Â 2D in case of rotation
                ALLOCATE( netcdf_data_2d(0:nx,1:ns) )
                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
                
                DO  k = 0, 2
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  Scalar grid points
                   IF ( k == 0 )  THEN
                      shift_x = 0.5 ; shift_y = 0.5
!                  
!--Â  Â  Â  Â  Â  Â  Â  Â  u grid points
                   ELSEIF ( k == 1 )  THEN
                      shift_x = 0.0 ; shift_y = 0.5
!                  
!--Â  Â  Â  Â  Â  Â  Â  Â  v grid points
                   ELSEIF ( k == 2 )  THEN
                      shift_x = 0.5 ; shift_y = 0.0
                   ENDIF

                   DO  j = 1, ns
                      IF( section(j,2) == -1 )  THEN
                         netcdf_data_2d(:,j) = -1.0_wp  ! section averaged along y
                      ELSE
                         DO  i = 0, nx
                            netcdf_data_2d(i,j) = init_model%origin_x          &
                                    + cos_ra * ( i + shift_x ) * dx            &
                                    + sin_ra * ( section(j,2) + shift_y ) * dy
                         ENDDO
                      ENDIF
                   ENDDO
                   
                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_eutm_xz(k,av),  &
                                           netcdf_data_2d, start = (/ 1, 1 /),   &
                                           count = (/ nx+1, ns /) )
                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
                   
                   DO  j = 1, ns
                      IF( section(j,2) == -1 )  THEN
                         netcdf_data_2d(:,j) = -1.0_wp  ! section averaged along y
                      ELSE
                         DO  i = 0, nx
                            netcdf_data_2d(i,j) = init_model%origin_y          &
                                    - sin_ra * ( i + shift_x ) * dx            &
                                    + cos_ra * ( section(j,2) + shift_y ) * dy
                         ENDDO
                      ENDIF
                   ENDDO
                   
                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_nutm_xz(k,av),  &
                                           netcdf_data_2d, start = (/ 1, 1 /),   &
                                           count = (/ nx+1, ns /) )
                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
                
                ENDDO
                DEALLOCATE( netcdf_data_2d )
             ENDIF
!
!--Â  Â  Â  Â  Â  Write lon and lat data
             ALLOCATE( lat(0:nx,1:ns) )
             ALLOCATE( lon(0:nx,1:ns) )
             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
             sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )

             DO  k = 0, 2
!               
!--Â  Â  Â  Â  Â  Â  Â Scalar grid points
                IF ( k == 0 )  THEN
                   shift_x = 0.5 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â u grid points
                ELSEIF ( k == 1 )  THEN
                   shift_x = 0.0 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â v grid points
                ELSEIF ( k == 2 )  THEN
                   shift_x = 0.5 ; shift_y = 0.0
                ENDIF

                DO  j = 1, ns
                   IF( section(j,2) == -1 )  THEN
                      lat(:,j) = -90.0_wp  ! section averaged along y
                      lon(:,j) = -180.0_wp  ! section averaged along y
                   ELSE
                      DO  i = 0, nx
                         eutm = init_model%origin_x            &
                              + cos_ra * ( i + shift_x ) * dx  &
                              + sin_ra * ( section(j,2) + shift_y ) * dy
                         nutm = init_model%origin_y            &
                              - sin_ra * ( i + shift_x ) * dx  &
                              + cos_ra * ( section(j,2) + shift_y ) * dy

                         CALL  convert_utm_to_geographic( crs_list,          &
                                                          eutm, nutm,        &
                                                          lon(i,j), lat(i,j) )
                      ENDDO
                   ENDIF
                ENDDO

                nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_lon_xz(k,av), &
                                     lon, start = (/ 1, 1 /),       &
                                     count = (/ nx+1, ns /) )
                CALL netcdf_handle_error( 'netcdf_define_header', 556 )

                nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_lat_xz(k,av), &
                                     lat, start = (/ 1, 1 /),       &
                                     count = (/ nx+1, ns /) )
                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
             ENDDO

             DEALLOCATE( lat )
             DEALLOCATE( lon )

          ENDIF


       CASE ( 'xz_ext' )

!
!--Â  Â  Â  Â 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_xz(av), NF90_GLOBAL, 'VAR_LIST', &
                                  var_list_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 168 )

          var_list = ';'
          i = 1
          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
             IF ( INDEX( do2d(av,i), 'xz' ) /= 0 )  THEN
                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
             ENDIF
             i = i + 1
          ENDDO

          IF ( av == 0 )  THEN
             var = '(xz)'
          ELSE
             var = '(xz_av)'
          ENDIF

          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
             message_string = 'netCDF file for cross-sections ' //           &
                              TRIM( var ) // ' from previous run found,' //  &
                              '&but this file cannot be extended due to' //  &
                              ' variable mismatch.' //                       &
                              '&New file is created instead.'
             CALL message( 'define_netcdf_header', 'PA0249', 0, 1, 0, 6, 0 )
             extend = .FALSE.
             RETURN
          ENDIF

!
!--Â  Â  Â  Â Calculate the number of current sections
          ns = 1
          DO WHILE ( section(ns,2) /= -9999  .AND.  ns <= 100 )
             ns = ns + 1
          ENDDO
          ns = ns - 1

!
!--Â  Â  Â  Â Get and compare the number of vertical cross sections
          nc_stat = NF90_INQ_VARID( id_set_xz(av), 'y_xz', id_var_y_xz(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 169 )

          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_y_xz(av), &
                                           dimids = id_dim_y_xz_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 170 )
          id_dim_y_xz(av) = id_dim_y_xz_old(1)

          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_y_xz(av), &
                                            len = ns_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 171 )

          IF ( ns /= ns_old )  THEN
             message_string = 'netCDF file for cross-sections ' //          &
                              TRIM( var ) // ' from previous run found,' // &
                              '&but this file cannot be extended due to' // &
                              ' mismatch in number of' //                   & 
                              ' cross sections.' //                         &
                              '&New file is created instead.'
             CALL message( 'define_netcdf_header', 'PA0250', 0, 1, 0, 6, 0 )
             extend = .FALSE.
             RETURN
          ENDIF

!
!--Â  Â  Â  Â Get and compare the heights of the cross sections
          ALLOCATE( netcdf_data(1:ns_old) )

          nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_y_xz(av), netcdf_data )
          CALL netcdf_handle_error( 'netcdf_define_header', 172 )

          DO  i = 1, ns
             IF ( section(i,2) /= -1 )  THEN
                IF ( ( ( section(i,2) + 0.5 ) * dy ) /= netcdf_data(i) )  THEN
                   message_string = 'netCDF file for cross-sections ' //       &
                               TRIM( var ) // ' from previous run found,' //   &
                               ' but this file cannot be extended' //          &
                               ' due to mismatch in cross' //                  &
                               ' section levels.' //                           &
                               ' New file is created instead.'
                   CALL message( 'define_netcdf_header', 'PA0251',             &
                                                                 0, 1, 0, 6, 0 )
                   extend = .FALSE.
                   RETURN
                ENDIF
             ELSE
                IF ( -1.0_wp /= netcdf_data(i) )  THEN
                   message_string = 'netCDF file for cross-sections ' //       &
                               TRIM( var ) // ' from previous run found,' //   &
                               ' but this file cannot be extended' //          &
                               ' due to mismatch in cross' //                  &
                               ' section levels.' //                           &
                               ' New file is created instead.'
                   CALL message( 'define_netcdf_header', 'PA0251',             &
                                                                 0, 1, 0, 6, 0 )
                   extend = .FALSE.
                   RETURN
                ENDIF
             ENDIF
          ENDDO

          DEALLOCATE( netcdf_data )

!
!--Â  Â  Â  Â Get the id of the time coordinate (unlimited coordinate) and its
!--Â  Â  Â  Â last index on the file. The next time level is do2d..count+1.
!--Â  Â  Â  Â The current time must be larger than the last output time
!--Â  Â  Â  Â on the file.
          nc_stat = NF90_INQ_VARID( id_set_xz(av), 'time', id_var_time_xz(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 173 )

          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_time_xz(av), &
                                           dimids = id_dim_time_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 174 )
          id_dim_time_xz(av) = id_dim_time_old(1)

          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_time_xz(av), &
                                            len = ntime_count )
          CALL netcdf_handle_error( 'netcdf_define_header', 175 )

!
!--Â  Â  Â  Â For non-parallel output use the last output time level of the netcdf
!--Â  Â  Â  Â file because the time dimension is unlimited. In case of parallel
!--Â  Â  Â  Â output the variable ntime_count could get the value of 9*10E36 because
!--Â  Â  Â  Â the time dimension is limited.
          IF ( netcdf_data_format < 5 ) do2d_xz_time_count(av) = ntime_count

          nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_time_xz(av),           &
                                  last_time_coordinate,                        &
                                  start = (/ do2d_xz_time_count(av) /),        &
                                  count = (/ 1 /) )
          CALL netcdf_handle_error( 'netcdf_define_header', 176 )

          IF ( last_time_coordinate(1) >= simulated_time )  THEN
             message_string = 'netCDF file for cross sections ' //             &
                              TRIM( var ) // ' from previous run found,' //    &
                              '&but this file cannot be extended becaus' //    &
                              'e the current output time' //                   &
                              '&is less or equal than the last output t' //    &
                              'ime on this file.' //                           &
                              '&New file is created instead.'
             CALL message( 'define_netcdf_header', 'PA0252', 0, 1, 0, 6, 0 )
             do2d_xz_time_count(av) = 0
             extend = .FALSE.
             RETURN
          ENDIF

          IF ( netcdf_data_format > 4 )  THEN
!
!--Â  Â  Â  Â  Â  Check if the needed number of output time levels is increased
!--Â  Â  Â  Â  Â  compared to the number of time levels in the existing file.
             IF ( ntdim_2d_xz(av) > ntime_count )  THEN
                message_string = 'netCDF file for cross sections ' // &
                                 TRIM( var ) // ' from previous run found,' // &
                                 '&but this file cannot be extended becaus' // &
                                 'e the number of output time levels has b' // &
                                 'een increased compared to the previous s' // &
                                 'imulation.' //                               &
                                 '&New file is created instead.'
                CALL message( 'define_netcdf_header', 'PA0390', 0, 1, 0, 6, 0 )
                do2d_xz_time_count(av) = 0
                extend = .FALSE.
!
!--Â  Â  Â  Â  Â  Â  Â Recalculate the needed time levels for the new file.
                IF ( av == 0 )  THEN
                   ntdim_2d_xz(0) = CEILING(                            &
                           ( end_time - MAX( skip_time_do2d_xz,         &
                                             simulated_time_at_begin )  &
                           ) / dt_do2d_xz )
                   IF ( do2d_at_begin )  ntdim_2d_xz(0) = ntdim_2d_xz(0) + 1
                ELSE
                   ntdim_2d_xz(1) = CEILING(                            &
                           ( end_time - MAX( skip_time_data_output_av,  &
                                             simulated_time_at_begin )  &
                           ) / dt_data_output_av )
                ENDIF
                RETURN
             ENDIF
          ENDIF

!
!--Â  Â  Â  Â Dataset seems to be extendable.
!--Â  Â  Â  Â Now get the variable ids.
          i = 1
          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
             IF ( INDEX( do2d(av,i), 'xz' ) /= 0 )  THEN
                nc_stat = NF90_INQ_VARID( id_set_xz(av), do2d(av,i), &
                                          id_var_do2d(av,i) )
                CALL netcdf_handle_error( 'netcdf_define_header', 177 )
#if defined( __netcdf4_parallel )
!
!--Â  Â  Â  Â  Â  Â  Â Set independent io operations for parallel io. Collective io
!--Â  Â  Â  Â  Â  Â  Â is only allowed in case of a 1d-decomposition along x, because
!--Â  Â  Â  Â  Â  Â  Â otherwise, not all PEs have output data.
                IF ( netcdf_data_format > 4 )  THEN
                   IF ( npey == 1 )  THEN
                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
                                                     id_var_do2d(av,i), &
                                                     NF90_COLLECTIVE )
                   ELSE
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Test simulations showed that the output of cross sections
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â by all PEs in data_output_2d using NF90_COLLECTIVE is 
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â faster than the output by the first row of PEs in
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â x-direction using NF90_INDEPENDENT.
                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
                                                     id_var_do2d(av,i), &
                                                     NF90_COLLECTIVE )
!                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
!                                                     id_var_do2d(av,i), &
!                                                     NF90_INDEPENDENT )
                   ENDIF
                   CALL netcdf_handle_error( 'netcdf_define_header', 455 )
                ENDIF
#endif
             ENDIF
             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_xz(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 433 )
          nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'title',         &
                                  TRIM( run_description_header ) //            &
                                  TRIM( time_average_text ) )
          CALL netcdf_handle_error( 'netcdf_define_header', 178 )
          nc_stat = NF90_ENDDEF( id_set_xz(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 434 )
          message_string = 'netCDF file for cross-sections ' //                &
                            TRIM( var ) // ' from previous run found.' //      &
                           '&This file will be extended.'
          CALL message( 'define_netcdf_header', 'PA0253', 0, 0, 0, 6, 0 )


       CASE ( 'yz_new' )

!
!--Â  Â  Â  Â Define some global attributes of the dataset
          IF ( av == 0 )  THEN
             CALL netcdf_create_global_atts( id_set_yz(av), 'yz', TRIM( run_description_header ), 179 )
             time_average_text = ' '
          ELSE
             CALL netcdf_create_global_atts( id_set_yz(av), 'yz_av', TRIM( run_description_header ), 179 )
             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
             nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'time_avg',   &
                                     TRIM( time_average_text ) )
             CALL netcdf_handle_error( 'netcdf_define_header', 180 )
          ENDIF

!
!--Â  Â  Â  Â Define time coordinate for yz sections.
!--Â  Â  Â  Â For parallel output the time dimensions has to be limited, otherwise
!--Â  Â  Â  Â the performance drops significantly.
          IF ( netcdf_data_format < 5 )  THEN
             CALL netcdf_create_dim( id_set_yz(av), 'time', NF90_UNLIMITED,    &
                                     id_dim_time_yz(av), 181 )
          ELSE
             CALL netcdf_create_dim( id_set_yz(av), 'time', ntdim_2d_yz(av),   &
                                     id_dim_time_yz(av), 526 )
          ENDIF

          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_time_yz(av) /),     &
                                  'time', NF90_DOUBLE, id_var_time_yz(av),     &
                                  'seconds since '//TRIM(init_model%origin_time), 'time', 182, 183, 000 )
          CALL netcdf_create_att( id_set_yz(av), id_var_time_yz(av), 'standard_name', 'time', 000)
          CALL netcdf_create_att( id_set_yz(av), id_var_time_yz(av), 'axis', 'T', 000)
!
!--Â  Â  Â  Â Define the spatial dimensions and coordinates for yz-sections.
!--Â  Â  Â  Â First, determine the number of vertical sections to be written.
          IF ( section(1,3) == -9999 )  THEN
             RETURN
          ELSE
             ns = 1
             DO WHILE ( section(ns,3) /= -9999  .AND.  ns <= 100 )
                ns = ns + 1
             ENDDO
             ns = ns - 1
          ENDIF

!
!--Â  Â  Â  Â Define x axis (for scalar position)
          CALL netcdf_create_dim( id_set_yz(av), 'x_yz', ns, id_dim_x_yz(av),  &
                                  184 )
          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_x_yz(av) /),        &
                                  'x_yz', NF90_DOUBLE, id_var_x_yz(av),        &
                                  'meters', '', 185, 186, 000 )
!
!--Â  Â  Â  Â Define x axis (for u position)
          CALL netcdf_create_dim( id_set_yz(av), 'xu_yz', ns,                  &
                                  id_dim_xu_yz(av), 377 )
          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_xu_yz(av) /),       &
                                  'xu_yz', NF90_DOUBLE, id_var_xu_yz(av),      &
                                  'meters', '', 378, 379, 000 )
!
!--Â  Â  Â  Â Define a variable to store the layer indices of the vertical cross
!--Â  Â  Â  Â sections
          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_x_yz(av) /),        &
                                  'ind_x_yz', NF90_DOUBLE,                     &
                                  id_var_ind_x_yz(av), 'gridpoints', '', 187,  &
                                  188, 000 )
!
!--Â  Â  Â  Â Define y-axis (for scalar position)
          CALL netcdf_create_dim( id_set_yz(av), 'y', ny+1, id_dim_y_yz(av),   &
                                  189 )
          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_y_yz(av) /), 'y',   &
                                  NF90_DOUBLE, id_var_y_yz(av), 'meters', '',  &
                                  190, 191, 000 )
!
!--Â  Â  Â  Â Define y-axis (for v position)
          CALL netcdf_create_dim( id_set_yz(av), 'yv', ny+1, id_dim_yv_yz(av), &
                                  380 )
          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_yv_yz(av) /), 'yv', &
                                  NF90_DOUBLE, id_var_yv_yz(av), 'meters', '', &
                                  381, 382, 000 )
!
!--Â  Â  Â  Â Define the two z-axes (zu and zw)
          CALL netcdf_create_dim( id_set_yz(av), 'zu', nz+2, id_dim_zu_yz(av), &
                                  192 )
          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_zu_yz(av) /), 'zu', &
                                  NF90_DOUBLE, id_var_zu_yz(av), 'meters', '', &
                                  193, 194, 000 )

          CALL netcdf_create_dim( id_set_yz(av), 'zw', nz+2, id_dim_zw_yz(av), &
                                  195 )
          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_zw_yz(av) /), 'zw', &
                                  NF90_DOUBLE, id_var_zw_yz(av), 'meters', '', &
                                  196, 197, 000 )
!
!--Â  Â  Â  Â Define UTM and geographic coordinates
          CALL define_geo_coordinates( id_set_yz(av),         &
                  (/ id_dim_x_yz(av), id_dim_xu_yz(av) /),    &
                  (/ id_dim_y_yz(av), id_dim_yv_yz(av) /),    &
                  id_var_eutm_yz(:,av), id_var_nutm_yz(:,av), &
                  id_var_lat_yz(:,av), id_var_lon_yz(:,av)    )
!
!--Â  Â  Â  Â Define coordinate-reference system
          CALL netcdf_create_crs( id_set_yz(av), 000 )

          IF ( land_surface )  THEN

             CALL netcdf_create_dim( id_set_yz(av), 'zs', nzs,                 &
                                     id_dim_zs_yz(av), 545 )
             CALL netcdf_create_var( id_set_yz(av), (/ id_dim_zs_yz(av) /),    &
                                     'zs', NF90_DOUBLE, id_var_zs_yz(av),      &
                                     'meters', '', 546, 547, 000 )

          ENDIF

!
!--Â  Â  Â  Â Define the variables
          var_list = ';'
          i = 1

          DO WHILE ( do2d(av,i)(1:1) /= ' ' )

             IF ( INDEX( do2d(av,i), 'yz' ) /= 0 )  THEN

!
!--Â  Â  Â  Â  Â  Â  Â Check for the grid
                found = .FALSE.
                SELECT CASE ( do2d(av,i) )
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Most variables are defined on the zu grid
                   CASE ( 'e_yz', 'nc_yz', 'nr_yz', 'p_yz', 'pc_yz',           &
                          'pr_yz','prr_yz', 'q_yz', 'qc_yz', 'ql_yz',          &
                          'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz',  &
                          's_yz',                                              &
                          'theta_yz', 'thetal_yz', 'thetav_yz' )

                      grid_x = 'x'
                      grid_y = 'y'
                      grid_z = 'zu'
!
!--Â  Â  Â  Â  Â  Â  Â  Â  u grid
                   CASE ( 'u_yz' )

                      grid_x = 'xu'
                      grid_y = 'y'
                      grid_z = 'zu'
!
!--Â  Â  Â  Â  Â  Â  Â  Â  v grid
                   CASE ( 'v_yz' )

                      grid_x = 'x'
                      grid_y = 'yv'
                      grid_z = 'zu'
!
!--Â  Â  Â  Â  Â  Â  Â  Â  w grid
                   CASE ( 'w_yz' )

                      grid_x = 'x'
                      grid_y = 'y'
                      grid_z = 'zw'


                   CASE DEFAULT
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for land surface quantities
                      IF ( land_surface )  THEN
                         CALL lsm_define_netcdf_grid( do2d(av,i), found,       &
                                                      grid_x, grid_y, grid_z )
                      ENDIF

                      IF ( .NOT. found )  THEN
                         CALL tcm_define_netcdf_grid( do2d(av,i), found,       &
                                                      grid_x, grid_y, grid_z )
                      ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for ocean quantities
                      IF ( .NOT. found  .AND.  ocean_mode )  THEN
                         CALL ocean_define_netcdf_grid( do2d(av,i), found,     &
                                                       grid_x, grid_y, grid_z )
                      ENDIF
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for radiation quantities
                      IF ( .NOT. found  .AND.  radiation )  THEN
                         CALL radiation_define_netcdf_grid( do2d(av,i), found, &
                                                            grid_x, grid_y,    &
                                                            grid_z )
                      ENDIF
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for SALSA quantities
                      IF ( .NOT. found  .AND.  salsa )  THEN
                         CALL salsa_define_netcdf_grid( do2d(av,i), found,     &
                                                        grid_x, grid_y, grid_z )
                      ENDIF                           
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for gust module quantities
                      IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
                         CALL gust_define_netcdf_grid( do2d(av,i), found,      &
                                                       grid_x, grid_y, grid_z )
                      ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for chemistry quantities
                      IF ( .NOT. found  .AND.  air_chemistry )  THEN
                         CALL chem_define_netcdf_grid( do2d(av,i), found,      &
                                                       grid_x, grid_y,         &
                                                       grid_z )
                      ENDIF
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Check for user-defined quantities
                      IF ( .NOT. found  .AND.  user_module_enabled )  THEN
                         CALL user_define_netcdf_grid( do2d(av,i), found,      &
                                                       grid_x, grid_y, grid_z )
                      ENDIF

                      IF ( .NOT. found )  THEN
                         WRITE ( message_string, * ) 'no grid defined for',    &
                                                ' variable ', TRIM( do2d(av,i) )
                         CALL message( 'define_netcdf_header', 'PA0244',       &
                                       0, 1, 0, 6, 0 )
                      ENDIF

                END SELECT

!
!--Â  Â  Â  Â  Â  Â  Â Select the respective dimension ids
                IF ( grid_x == 'x' )  THEN
                   id_x = id_dim_x_yz(av)
                ELSEIF ( grid_x == 'xu' )  THEN
                   id_x = id_dim_xu_yz(av)
                ENDIF

                IF ( grid_y == 'y' )  THEN
                   id_y = id_dim_y_yz(av)
                ELSEIF ( grid_y == 'yv' )  THEN
                   id_y = id_dim_yv_yz(av)
                ENDIF

                IF ( grid_z == 'zu' )  THEN
                   id_z = id_dim_zu_yz(av)
                ELSEIF ( grid_z == 'zw' )  THEN
                   id_z = id_dim_zw_yz(av)
                ELSEIF ( grid_z == 'zs' )  THEN
                   id_z = id_dim_zs_yz(av)
                ENDIF

!
!--Â  Â  Â  Â  Â  Â  Â Define the grid
                CALL netcdf_create_var( id_set_yz(av),  (/ id_x, id_y, id_z,   &
                                        id_dim_time_yz(av) /), do2d(av,i),     &
                                        nc_precision(3), id_var_do2d(av,i),    &
                                        TRIM( do2d_unit(av,i) ), do2d(av,i),   &
                                        198, 199, 356, .TRUE. )

#if defined( __netcdf4_parallel )
                IF ( netcdf_data_format > 4 )  THEN
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Set no fill for every variable to increase performance.
                   nc_stat = NF90_DEF_VAR_FILL( id_set_yz(av),     &
                                                id_var_do2d(av,i), &
                                                1, 0 )
                   CALL netcdf_handle_error( 'netcdf_define_header', 535 )
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Set independent io operations for parallel io. Collective io
!--Â  Â  Â  Â  Â  Â  Â  Â  is only allowed in case of a 1d-decomposition along y,
!--Â  Â  Â  Â  Â  Â  Â  Â  because otherwise, not all PEs have output data.
                   IF ( npex == 1 )  THEN
                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av),     &
                                                     id_var_do2d(av,i), &
                                                     NF90_COLLECTIVE )
                   ELSE
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Test simulations showed that the output of cross sections
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â by all PEs in data_output_2d using NF90_COLLECTIVE is 
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â faster than the output by the first row of PEs in
!--Â  Â  Â  Â  Â  Â  Â  Â  Â  Â y-direction using NF90_INDEPENDENT.
                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av),     &
                                                     id_var_do2d(av,i), &
                                                     NF90_COLLECTIVE )
!                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av),     &
!                                                     id_var_do2d(av,i), &
!                                                     NF90_INDEPENDENT )
                   ENDIF
                   CALL netcdf_handle_error( 'netcdf_define_header', 450 )
                ENDIF
#endif
                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'

             ENDIF

             i = i + 1

          ENDDO

!
!--Â  Â  Â  Â No arrays to output. Close the netcdf file and return.
          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_yz(av), NF90_GLOBAL, 'VAR_LIST', &
                                  var_list )
          CALL netcdf_handle_error( 'netcdf_define_header', 200 )

!
!--Â  Â  Â  Â Set general no fill, otherwise the performance drops significantly for
!--Â  Â  Â  Â parallel output.
          nc_stat = NF90_SET_FILL( id_set_yz(av), NF90_NOFILL, oldmode )
          CALL netcdf_handle_error( 'netcdf_define_header', 531 )

!
!--Â  Â  Â  Â Leave netCDF define mode
          nc_stat = NF90_ENDDEF( id_set_yz(av) )
          CALL netcdf_handle_error( 'netcdf_define_header', 201 )

!
!--Â  Â  Â  Â These data are only written by PE0 for parallel output to increase
!--Â  Â  Â  Â the performance.
          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN

!
!--Â  Â  Â  Â  Â  Write axis data: x_yz, y, zu, zw
             ALLOCATE( netcdf_data(1:ns) )

!
!--Â  Â  Â  Â  Â  Write x_yz data (shifted by +dx/2)
             DO  i = 1, ns
                IF( section(i,3) == -1 )  THEN
                   netcdf_data(i) = -1.0_wp  ! section averaged along x
                ELSE
                   netcdf_data(i) = ( section(i,3) + 0.5_wp ) * dx
                ENDIF
             ENDDO
             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_x_yz(av), &
                                     netcdf_data, start = (/ 1 /),   &
                                     count = (/ ns /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 202 )

!
!--Â  Â  Â  Â  Â  Write x_yz data (xu grid)
             DO  i = 1, ns
                IF( section(i,3) == -1 )  THEN
                   netcdf_data(i) = -1.0_wp  ! section averaged along x
                ELSE
                   netcdf_data(i) = section(i,3) * dx
                ENDIF
             ENDDO
             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_xu_yz(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ ns /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 383 )

!
!--Â  Â  Â  Â  Â  Write gridpoint number data
             netcdf_data(1:ns) = section(1:ns,3)
             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_ind_x_yz(av), &
                                     netcdf_data, start = (/ 1 /),       &
                                     count = (/ ns /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 203 )

             DEALLOCATE( netcdf_data )

!
!--Â  Â  Â  Â  Â  Write data for y (shifted by +dy/2) and yv axis
             ALLOCATE( netcdf_data(0:ny) )

             DO  j = 0, ny
                netcdf_data(j) = ( j + 0.5_wp ) * dy
             ENDDO

             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_y_yz(av), &
                                     netcdf_data, start = (/ 1 /),   &
                                     count = (/ ny+1 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 204 )

             DO  j = 0, ny
                netcdf_data(j) = j * dy
             ENDDO

             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_yv_yz(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ ny+1 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 384 )

             DEALLOCATE( netcdf_data )

!
!--Â  Â  Â  Â  Â  Write zu and zw data (vertical axes)
             ALLOCATE( netcdf_data(0:nz+1) )

             netcdf_data(0:nz+1) = zu(nzb:nzt+1)
             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zu_yz(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ nz+2 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 205 )

             netcdf_data(0:nz+1) = zw(nzb:nzt+1)
             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zw_yz(av), &
                                     netcdf_data, start = (/ 1 /),    &
                                     count = (/ nz+2 /) )
             CALL netcdf_handle_error( 'netcdf_define_header', 206 )

             DEALLOCATE( netcdf_data )
!
!--Â  Â  Â  Â  Â  Write UTM coordinates
             IF ( init_model%rotation_angle == 0.0_wp )  THEN
!
!--Â  Â  Â  Â  Â  Â  Â 1D in case of no rotation
                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
!
!--Â  Â  Â  Â  Â  Â  Â x coordinates
                ALLOCATE( netcdf_data(1:ns) )
                DO  k = 0, 2
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  Scalar grid points
                   IF ( k == 0 )  THEN
                      shift_x = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  u grid points
                   ELSEIF ( k == 1 )  THEN
                      shift_x = 0.0
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  v grid points
                   ELSEIF ( k == 2 )  THEN
                      shift_x = 0.5
                   ENDIF
                
                   DO  i = 1, ns
                      IF( section(i,3) == -1 )  THEN
                         netcdf_data(i) = -1.0_wp  ! section averaged along x
                      ELSE
                         netcdf_data(i) = init_model%origin_x &
                                     + cos_ra * ( section(i,3) + shift_x ) * dx
                      ENDIF
                   ENDDO
                
                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_eutm_yz(k,av),&
                                           netcdf_data, start = (/ 1 /),   &
                                           count = (/ ns /) )
                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )

                ENDDO
                DEALLOCATE( netcdf_data )
!
!--Â  Â  Â  Â  Â  Â  Â y coordinates
                ALLOCATE( netcdf_data(0:ny) )
                DO  k = 0, 2
!
!--Â  Â  Â  Â  Â  Â  Â  Â  Scalar grid points
                   IF ( k == 0 )  THEN
                      shift_y = 0.5
!
!--Â  Â  Â  Â  Â  Â  Â  Â  u grid points
                   ELSEIF ( k == 1 )  THEN
                      shift_y = 0.5
!
!--Â  Â  Â  Â  Â  Â  Â  Â  v grid points
                   ELSEIF ( k == 2 )  THEN
                      shift_y = 0.0
                   ENDIF

                   DO  i = 0, ny
                     netcdf_data(i) = init_model%origin_y            &
                                    + cos_ra * ( i + shift_y ) * dy
                   ENDDO

                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_nutm_yz(k,av),&
                                           netcdf_data, start = (/ 1 /),   &
                                           count = (/ ny+1 /) )
                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )

                ENDDO
                DEALLOCATE( netcdf_data )

             ELSE
!
!--Â  Â  Â  Â  Â  Â  Â 2D in case of rotation
                ALLOCATE( netcdf_data_2d(1:ns,0:ny) )
                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )

                DO  k = 0, 2
!               
!--Â  Â  Â  Â  Â  Â  Â  Â  Scalar grid points
                   IF ( k == 0 )  THEN
                      shift_x = 0.5 ; shift_y = 0.5
!                  
!--Â  Â  Â  Â  Â  Â  Â  Â  u grid points
                   ELSEIF ( k == 1 )  THEN
                      shift_x = 0.0 ; shift_y = 0.5
!                  
!--Â  Â  Â  Â  Â  Â  Â  Â  v grid points
                   ELSEIF ( k == 2 )  THEN
                      shift_x = 0.5 ; shift_y = 0.0
                   ENDIF

                   DO  j = 0, ny
                      DO  i = 1, ns
                         IF( section(i,3) == -1 )  THEN
                            netcdf_data_2d(i,:) = -1.0_wp !section averaged along x
                         ELSE
                            netcdf_data_2d(i,j) = init_model%origin_x          &
                                    + cos_ra * ( section(i,3) + shift_x ) * dx &
                                    + sin_ra * ( j + shift_y ) * dy
                         ENDIF
                      ENDDO
                   ENDDO
                   
                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_eutm_yz(k,av),  &
                                           netcdf_data_2d, start = (/ 1, 1 /),   &
                                           count = (/ ns, ny+1 /) )
                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
                   
                   DO  j = 0, ny
                      DO  i = 1, ns
                         IF( section(i,3) == -1 )  THEN
                            netcdf_data_2d(i,:) = -1.0_wp !section averaged along x
                         ELSE
                            netcdf_data_2d(i,j) = init_model%origin_y          &
                                    - sin_ra * ( section(i,3) + shift_x ) * dx &
                                    + cos_ra * ( j + shift_y ) * dy
                         ENDIF
                      ENDDO
                   ENDDO

                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_nutm_yz(k,av),  &
                                           netcdf_data_2d, start = (/ 1, 1 /),   &
                                           count = (/ ns, ny+1 /) )
                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
                
                ENDDO
                DEALLOCATE( netcdf_data_2d )
             ENDIF
!
!--Â  Â  Â  Â  Â  Write lon and lat data
             ALLOCATE( lat(1:ns,0:ny) )
             ALLOCATE( lon(1:ns,0:ny) )
             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
             sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )

             DO  k = 0, 2
!               
!--Â  Â  Â  Â  Â  Â  Â Scalar grid points
                IF ( k == 0 )  THEN
                   shift_x = 0.5 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â u grid points
                ELSEIF ( k == 1 )  THEN
                   shift_x = 0.0 ; shift_y = 0.5
!               
!--Â  Â  Â  Â  Â  Â  Â v grid points
                ELSEIF ( k == 2 )  THEN
                   shift_x = 0.5 ; shift_y = 0.0
                ENDIF

                DO  j = 0, ny
                   DO  i = 1, ns
                      IF( section(i,3) == -1 )  THEN
                         lat(i,:) = -90.0_wp   ! section averaged along x
                         lon(i,:) = -180.0_wp  ! section averaged along x
                      ELSE
                         eutm = init_model%origin_x            &
                              + cos_ra * ( section(i,3) + shift_x ) * dx  &
                              + sin_ra * ( j + shift_y ) * dy
                         nutm = init_model%origin_y            &
                              - sin_ra * ( section(i,3) + shift_x ) * dx  &
                              + cos_ra * ( j + shift_y ) * dy

                         CALL  convert_utm_to_geographic( crs_list,          &
                                                          eutm, nutm,        &
                                                          lon(i,j), lat(i,j) )
                      ENDIF
                   ENDDO
                ENDDO

                nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_lon_yz(k,av), &
                                     lon, start = (/ 1, 1 /),       &
                                     count = (/ ns, ny+1 /) )
                CALL netcdf_handle_error( 'netcdf_define_header', 556 )

                nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_lat_yz(k,av), &
                                     lat, start = (/ 1, 1 /),       &
                                     count = (/ ns, ny+1 /) )
                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
             ENDDO

             DEALLOCATE( lat )
             DEALLOCATE( lon )

          ENDIF


       CASE ( 'yz_ext' )

!
!--Â  Â  Â  Â 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_yz(av), NF90_GLOBAL, 'VAR_LIST', &
                                  var_list_old )
          CALL netcdf_handle_error( 'netcdf_define_header', 207 )

          var_list = ';'
          i = 1
          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
             IF ( INDEX( do2d(av,i), 'yz' ) /= 0 )  THEN
                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
             ENDIF
             i = i + 1
          ENDDO

Â