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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

  • Property svn:keywords set to Id
File size: 319.0 KB
Line 
1!> @file netcdf_interface_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: netcdf_interface_mod.f90 4180 2019-08-21 14:37:54Z scharf $
27! -Introduce new vertical dimension for plant-canopy output.
28! -Temporarlily disable masked output for soil (merge from branch resler)
29!
30! 4069 2019-07-01 14:05:51Z Giersch
31! Masked output running index mid has been introduced as a local variable to
32! avoid runtime error (Loop variable has been modified) in time_integration
33!
34! 4046 2019-06-21 17:32:04Z knoop
35! removal of special treatment for usm_define_netcdf_grid call
36!
37! 4039 2019-06-18 10:32:41Z suehring
38! Rename subroutines in module for diagnostic quantities
39!
40! 4029 2019-06-14 14:04:35Z raasch
41! netcdf variable NF90_NOFILL is used as argument instead of "1" in calls to NF90_DEF_VAR_FILL
42!
43! 3995 2019-05-22 18:59:54Z suehring
44! output of turbulence intensity added
45!
46! 3994 2019-05-22 18:08:09Z suehring
47! remove origin time from time unit, compose origin_time_string within
48! subroutine netcdf_create_global_atts
49!
50! 3954 2019-05-06 12:49:42Z gronemeier
51! bugfix: corrected format for date_time_string
52!
53! 3953 2019-05-06 12:11:55Z gronemeier
54! bugfix: set origin_time and starting point of time coordinate according to
55!         day_of_year_init and time_utc_init
56!
57! 3942 2019-04-30 13:08:30Z kanani
58! Add specifier to netcdf_handle_error to simplify identification of attribute
59! causing the error
60!
61! 3766 2019-02-26 16:23:41Z raasch
62! bugfix in im_define_netcdf_grid argument list
63!
64! 3745 2019-02-15 18:57:56Z suehring
65! Add indoor model
66!
67! 3744 2019-02-15 18:38:58Z suehring
68! Bugfix: - initialize return values to ensure they are set before returning
69!           (routine define_geo_coordinates)
70!         - change order of dimensions for some variables
71!
72! 3727 2019-02-08 14:52:10Z gronemeier
73! make several routines publicly available
74!
75! 3701 2019-01-26 18:57:21Z knoop
76! Statement added to prevent compiler warning about unused variable
77!
78! 3655 2019-01-07 16:51:22Z knoop
79! Move the control parameter "salsa" from salsa_mod to control_parameters
80! (M. Kurppa)
81!
82!
83! Description:
84! ------------
85!> In case of extend = .FALSE.:
86!> Define all necessary dimensions, axes and variables for the different
87!> netCDF datasets. This subroutine is called from check_open after a new
88!> dataset is created. It leaves the open netCDF files ready to write.
89!>
90!> In case of extend = .TRUE.:
91!> Find out if dimensions and variables of an existing file match the values
92!> of the actual run. If so, get all necessary information (ids, etc.) from
93!> this file.
94!>
95!> Parameter av can assume values 0 (non-averaged data) and 1 (time averaged
96!> data)
97!>
98!> @todo calculation of output time levels for parallel NetCDF still does not
99!>       cover every exception (change of dt_do, end_time in restart)
100!> @todo timeseries and profile output still needs to be rewritten to allow
101!>       modularization
102!> @todo output 2d UTM coordinates without global arrays
103!> @todo output longitude/latitude also with non-parallel output (3d and xy)
104!------------------------------------------------------------------------------!
105 MODULE netcdf_interface
106
107    USE control_parameters,                                                    &
108        ONLY:  biometeorology, fl_max,                                         &
109               max_masks, multi_agent_system_end,                              &
110               multi_agent_system_start, var_fl_max, varnamelength
111    USE kinds
112#if defined( __netcdf )
113    USE NETCDF
114#endif
115    USE mas_global_attributes,                                                 &
116        ONLY:  dim_size_agtnum
117
118    USE netcdf_data_input_mod,                                                 &
119        ONLY: coord_ref_sys, init_model
120
121    PRIVATE
122
123    CHARACTER (LEN=16), DIMENSION(13) ::  agt_var_names =                      &
124          (/ 'ag_id           ', 'ag_x            ', 'ag_y            ',       &
125             'ag_wind         ', 'ag_temp         ', 'ag_group        ',       &
126             'PM10            ', 'PM25            ', 'ag_iPT          ',       &
127             'ag_uv           ', 'not_used        ', 'not_used        ',       &
128             'not_used        ' /)
129
130    CHARACTER (LEN=16), DIMENSION(13) ::  agt_var_units = &
131          (/ 'dim_less        ', 'meters          ', 'meters          ',       &
132             'm/s             ', 'K               ', 'dim_less        ',       &
133             'tbd             ', 'tbd             ', 'tbd             ',       &
134             'tbd             ', 'not_used        ', 'not_used        ',       &
135             'not_used        ' /)
136
137    INTEGER(iwp), PARAMETER ::  dopr_norm_num = 7, dopts_num = 29, dots_max = 100
138
139    CHARACTER (LEN=7), DIMENSION(dopr_norm_num) ::  dopr_norm_names =          &
140         (/ 'wtheta0', 'ws2    ', 'tsw2   ', 'ws3    ', 'ws2tsw ', 'wstsw2 ',  &
141            'z_i    ' /)
142
143    CHARACTER (LEN=7), DIMENSION(dopr_norm_num) ::  dopr_norm_longnames =      &
144         (/ 'wtheta0', 'w*2    ', 't*w2   ', 'w*3    ', 'w*2t*w ', 'w*t*w2 ',  &
145            'z_i    ' /)
146
147    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_label =                   &
148          (/ 'tnpt   ', 'x_     ', 'y_     ', 'z_     ', 'z_abs  ', 'u      ', &
149             'v      ', 'w      ', 'u"     ', 'v"     ', 'w"     ', 'npt_up ', &
150             'w_up   ', 'w_down ', 'radius ', 'r_min  ', 'r_max  ', 'npt_max', &
151             'npt_min', 'x*2    ', 'y*2    ', 'z*2    ', 'u*2    ', 'v*2    ', &
152             'w*2    ', 'u"2    ', 'v"2    ', 'w"2    ', 'npt*2  ' /)
153
154    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_unit =                    &
155          (/ 'number ', 'm      ', 'm      ', 'm      ', 'm      ', 'm/s    ', &
156             'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'number ', &
157             'm/s    ', 'm/s    ', 'm      ', 'm      ', 'm      ', 'number ', &
158             'number ', 'm2     ', 'm2     ', 'm2     ', 'm2/s2  ', 'm2/s2  ', &
159             'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'number2' /)
160
161    INTEGER(iwp) ::  dots_num  = 25  !< number of timeseries defined by default
162    INTEGER(iwp) ::  dots_soil = 26  !< starting index for soil-timeseries
163    INTEGER(iwp) ::  dots_rad  = 32  !< starting index for radiation-timeseries
164
165    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label =                    &
166          (/ 'E            ', 'E*           ', 'dt           ',                &
167             'us*          ', 'th*          ', 'umax         ',                &
168             'vmax         ', 'wmax         ', 'div_new      ',                &
169             'div_old      ', 'zi_wtheta    ', 'zi_theta     ',                &
170             'w*           ', 'w"theta"0    ', 'w"theta"     ',                &
171             'wtheta       ', 'theta(0)     ', 'theta(z_mo)  ',                &
172             'w"u"0        ', 'w"v"0        ', 'w"q"0        ',                &
173             'ol           ', 'q*           ', 'w"s"         ',                &
174             's*           ', 'ghf          ', 'qsws_liq     ',                &
175             'qsws_soil    ', 'qsws_veg     ', 'r_a          ',                &
176             'r_s          ',                                                  &
177             'rad_net      ', 'rad_lw_in    ', 'rad_lw_out   ',                &
178             'rad_sw_in    ', 'rad_sw_out   ', 'rrtm_aldif   ',                &
179             'rrtm_aldir   ', 'rrtm_asdif   ', 'rrtm_asdir   ',                &
180             ( 'unknown      ', i9 = 1, dots_max-40 ) /)
181
182    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit =                     &
183          (/ 'm2/s2        ', 'm2/s2        ', 's            ',                &
184             'm/s          ', 'K            ', 'm/s          ',                &
185             'm/s          ', 'm/s          ', 's-1          ',                &
186             's-1          ', 'm            ', 'm            ',                &
187             'm/s          ', 'K m/s        ', 'K m/s        ',                &
188             'K m/s        ', 'K            ', 'K            ',                &
189             'm2/s2        ', 'm2/s2        ', 'kg m/s       ',                &
190             'm            ', 'kg/kg        ', 'kg m/(kg s)  ',                &
191             'kg/kg        ', 'W/m2         ', 'W/m2         ',                &
192             'W/m2         ', 'W/m2         ', 's/m          ',                &
193             's/m          ',                                                  &
194             'W/m2         ', 'W/m2         ', 'W/m2         ',                &
195             'W/m2         ', 'W/m2         ', '             ',                &
196             '             ', '             ', '             ',                &
197             ( 'unknown      ', i9 = 1, dots_max-40 ) /)
198
199    CHARACTER (LEN=16) :: heatflux_output_unit     !< unit for heatflux output
200    CHARACTER (LEN=16) :: waterflux_output_unit    !< unit for waterflux output
201    CHARACTER (LEN=16) :: momentumflux_output_unit !< unit for momentumflux output
202
203    CHARACTER (LEN=9), DIMENSION(300) ::  dopr_unit = 'unknown'
204
205    CHARACTER (LEN=7), DIMENSION(0:1,500) ::  do2d_unit, do3d_unit
206
207!    CHARACTER (LEN=16), DIMENSION(25) ::  prt_var_names = &
208!          (/ 'pt_age          ', 'pt_dvrp_size    ', 'pt_origin_x     ', &
209!             'pt_origin_y     ', 'pt_origin_z     ', 'pt_radius       ', &
210!             'pt_speed_x      ', 'pt_speed_y      ', 'pt_speed_z      ', &
211!             'pt_weight_factor', 'pt_x            ', 'pt_y            ', &
212!             'pt_z            ', 'pt_color        ', 'pt_group        ', &
213!             'pt_tailpoints   ', 'pt_tail_id      ', 'pt_density_ratio', &
214!             'pt_exp_arg      ', 'pt_exp_term     ', 'not_used        ', &
215!             'not_used        ', 'not_used        ', 'not_used        ', &
216!             'not_used        ' /)
217
218!    CHARACTER (LEN=16), DIMENSION(25) ::  prt_var_units = &
219!          (/ 'seconds         ', 'meters          ', 'meters          ', &
220!             'meters          ', 'meters          ', 'meters          ', &
221!             'm/s             ', 'm/s             ', 'm/s             ', &
222!             'factor          ', 'meters          ', 'meters          ', &
223!             'meters          ', 'none            ', 'none            ', &
224!             'none            ', 'none            ', 'ratio           ', &
225!             'none            ', 'none            ', 'not_used        ', &
226!             'not_used        ', 'not_used        ', 'not_used        ', &
227!             'not_used        ' /)
228
229    CHARACTER(LEN=20), DIMENSION(11) ::  netcdf_precision = ' '
230    CHARACTER(LEN=40) ::  netcdf_data_format_string
231
232    INTEGER(iwp) ::  id_dim_agtnum, id_dim_time_agt,                           &
233                     id_dim_time_fl, id_dim_time_pr,                           &
234                     id_dim_time_pts, id_dim_time_sp, id_dim_time_ts,          &
235                     id_dim_x_sp, id_dim_y_sp, id_dim_zu_sp, id_dim_zw_sp,     &
236                     id_set_agt, id_set_fl, id_set_pr, id_set_prt, id_set_pts, &
237                     id_set_sp, id_set_ts, id_var_agtnum, id_var_time_agt,     &
238                     id_var_time_fl, id_var_rnoa_agt, id_var_time_pr,          &
239                     id_var_time_pts, id_var_time_sp, id_var_time_ts,          &
240                     id_var_x_sp, id_var_y_sp, id_var_zu_sp, id_var_zw_sp,     &
241                     nc_stat
242
243
244    INTEGER(iwp), DIMENSION(0:1) ::  id_dim_time_xy, id_dim_time_xz, &
245                    id_dim_time_yz, id_dim_time_3d, id_dim_x_xy, id_dim_xu_xy, &
246                    id_dim_x_xz, id_dim_xu_xz, id_dim_x_yz, id_dim_xu_yz, &
247                    id_dim_x_3d, id_dim_xu_3d, id_dim_y_xy, id_dim_yv_xy, &
248                    id_dim_y_xz, id_dim_yv_xz, id_dim_y_yz, id_dim_yv_yz, &
249                    id_dim_y_3d, id_dim_yv_3d, id_dim_zs_xy, id_dim_zs_xz, &
250                    id_dim_zs_yz, id_dim_zs_3d, id_dim_zpc_3d, &
251                    id_dim_zu_xy, id_dim_zu1_xy, &
252                    id_dim_zu_xz, id_dim_zu_yz, id_dim_zu_3d, id_dim_zw_xy, &
253                    id_dim_zw_xz, id_dim_zw_yz, id_dim_zw_3d, id_set_xy, &
254                    id_set_xz, id_set_yz, id_set_3d, id_var_ind_x_yz, &
255                    id_var_ind_y_xz, id_var_ind_z_xy, id_var_time_xy, &
256                    id_var_time_xz, id_var_time_yz, id_var_time_3d, id_var_x_xy, &
257                    id_var_xu_xy, id_var_x_xz, id_var_xu_xz, id_var_x_yz, &
258                    id_var_xu_yz, id_var_x_3d, id_var_xu_3d, id_var_y_xy, &
259                    id_var_yv_xy, id_var_y_xz, id_var_yv_xz, id_var_y_yz, &
260                    id_var_yv_yz, id_var_y_3d, id_var_yv_3d, id_var_zs_xy, &
261                    id_var_zs_xz, id_var_zs_yz, id_var_zs_3d, id_var_zpc_3d, &
262                    id_var_zusi_xy, id_var_zusi_3d, id_var_zu_xy, id_var_zu1_xy, id_var_zu_xz, &
263                    id_var_zu_yz, id_var_zu_3d, id_var_zwwi_xy, id_var_zwwi_3d, &
264                    id_var_zw_xy, id_var_zw_xz, id_var_zw_yz, id_var_zw_3d
265
266    INTEGER(iwp), DIMENSION(0:2,0:1) ::  id_var_eutm_3d, id_var_nutm_3d, &
267                                         id_var_eutm_xy, id_var_nutm_xy, &
268                                         id_var_eutm_xz, id_var_nutm_xz, &
269                                         id_var_eutm_yz, id_var_nutm_yz
270
271    INTEGER(iwp), DIMENSION(0:2,0:1) ::  id_var_lat_3d, id_var_lon_3d, &
272                                         id_var_lat_xy, id_var_lon_xy, &
273                                         id_var_lat_xz, id_var_lon_xz, &
274                                         id_var_lat_yz, id_var_lon_yz
275
276    INTEGER ::  netcdf_data_format = 2  !< NetCDF3 64bit offset format
277    INTEGER ::  netcdf_deflate = 0      !< NetCDF compression, default: no
278                                        !< compression
279
280    INTEGER(iwp)                 ::  dofl_time_count
281    INTEGER(iwp), DIMENSION(10)  ::  id_var_dospx, id_var_dospy
282    INTEGER(iwp), DIMENSION(20)  ::  id_var_agt
283!    INTEGER(iwp), DIMENSION(20)  ::  id_var_prt
284    INTEGER(iwp), DIMENSION(11)  ::  nc_precision
285    INTEGER(iwp), DIMENSION(dopr_norm_num) ::  id_var_norm_dopr
286
287    INTEGER(iwp), DIMENSION(fl_max) ::  id_dim_x_fl, id_dim_y_fl, id_dim_z_fl
288    INTEGER(iwp), DIMENSION(fl_max) ::  id_var_x_fl, id_var_y_fl, id_var_z_fl
289
290    CHARACTER (LEN=20), DIMENSION(fl_max*var_fl_max) :: dofl_label
291    CHARACTER (LEN=20), DIMENSION(fl_max*var_fl_max) :: dofl_unit
292    CHARACTER (LEN=20), DIMENSION(fl_max) :: dofl_dim_label_x
293    CHARACTER (LEN=20), DIMENSION(fl_max) :: dofl_dim_label_y
294    CHARACTER (LEN=20), DIMENSION(fl_max) :: dofl_dim_label_z
295
296    INTEGER(iwp), DIMENSION(fl_max*var_fl_max) :: id_var_dofl
297
298    INTEGER(iwp), DIMENSION(dopts_num,0:10) ::  id_var_dopts
299    INTEGER(iwp), DIMENSION(0:1,500)        ::  id_var_do2d, id_var_do3d
300    INTEGER(iwp), DIMENSION(100,0:99)       ::  id_dim_z_pr, id_var_dopr, &
301                                                id_var_z_pr
302    INTEGER(iwp), DIMENSION(dots_max,0:99)  ::  id_var_dots
303
304!
305!-- Masked output
306    CHARACTER (LEN=7), DIMENSION(max_masks,0:1,100) ::  domask_unit
307
308    LOGICAL ::  output_for_t0 = .FALSE.
309
310    INTEGER(iwp), DIMENSION(1:max_masks,0:1) ::  id_dim_time_mask, id_dim_x_mask, &
311                   id_dim_xu_mask, id_dim_y_mask, id_dim_yv_mask, id_dim_zs_mask, &
312                   id_dim_zu_mask, id_dim_zw_mask, &
313                   id_set_mask, &
314                   id_var_time_mask, id_var_x_mask, id_var_xu_mask, &
315                   id_var_y_mask, id_var_yv_mask, id_var_zs_mask, &
316                   id_var_zu_mask, id_var_zw_mask, &
317                   id_var_zusi_mask, id_var_zwwi_mask
318
319    INTEGER(iwp), DIMENSION(0:2,1:max_masks,0:1) ::  id_var_eutm_mask, &
320                                                     id_var_nutm_mask
321
322    INTEGER(iwp), DIMENSION(0:2,1:max_masks,0:1) ::  id_var_lat_mask, &
323                                                     id_var_lon_mask
324
325    INTEGER(iwp), DIMENSION(1:max_masks,0:1,100) ::  id_var_domask
326
327    REAL(wp) ::  fill_value = -9999.0_wp    !< value for the _FillValue attribute
328
329
330    PUBLIC  dofl_dim_label_x, dofl_dim_label_y, dofl_dim_label_z, dofl_label,  &
331            dofl_time_count, dofl_unit, domask_unit, dopr_unit, dopts_num,     &
332            dots_label, dots_max, dots_num, dots_rad, dots_soil, dots_unit,    &
333            do2d_unit, do3d_unit, fill_value, id_set_agt, id_set_fl,           &
334            id_set_mask, id_set_pr, id_set_prt, id_set_pts, id_set_sp,         &
335            id_set_ts, id_set_xy, id_set_xz, id_set_yz, id_set_3d, id_var_agt, &
336            id_var_domask, id_var_dofl, id_var_dopr, id_var_dopts,             &
337            id_var_dospx, id_var_dospy, id_var_dots, id_var_do2d, id_var_do3d, &
338            id_var_norm_dopr, id_var_time_agt, id_var_time_fl,                 &
339            id_var_time_mask, id_var_time_pr, id_var_rnoa_agt, id_var_time_pts,&
340            id_var_time_sp, id_var_time_ts,                                    &
341            id_var_time_xy, id_var_time_xz, id_var_time_yz, id_var_time_3d,    &
342            id_var_x_fl, id_var_y_fl, id_var_z_fl,  nc_stat,                   &
343            netcdf_data_format, netcdf_data_format_string, netcdf_deflate,     &
344            netcdf_precision, output_for_t0, heatflux_output_unit,             &
345            waterflux_output_unit, momentumflux_output_unit
346
347    SAVE
348
349    INTERFACE netcdf_create_dim
350       MODULE PROCEDURE netcdf_create_dim
351    END INTERFACE netcdf_create_dim
352
353    INTERFACE netcdf_create_file
354       MODULE PROCEDURE netcdf_create_file
355    END INTERFACE netcdf_create_file
356
357    INTERFACE netcdf_create_global_atts
358       MODULE PROCEDURE netcdf_create_global_atts
359    END INTERFACE netcdf_create_global_atts
360
361    INTERFACE netcdf_create_var
362       MODULE PROCEDURE netcdf_create_var
363    END INTERFACE netcdf_create_var
364
365    INTERFACE netcdf_create_att
366       MODULE PROCEDURE netcdf_create_att_string
367    END INTERFACE netcdf_create_att
368
369    INTERFACE netcdf_define_header
370       MODULE PROCEDURE netcdf_define_header
371    END INTERFACE netcdf_define_header
372
373    INTERFACE netcdf_handle_error
374       MODULE PROCEDURE netcdf_handle_error
375    END INTERFACE netcdf_handle_error
376
377    INTERFACE netcdf_open_write_file
378       MODULE PROCEDURE netcdf_open_write_file
379    END INTERFACE netcdf_open_write_file
380
381    PUBLIC netcdf_create_att, netcdf_create_dim, netcdf_create_file,           &
382           netcdf_create_global_atts, netcdf_create_var, netcdf_define_header, &
383           netcdf_handle_error, netcdf_open_write_file
384
385 CONTAINS
386
387 SUBROUTINE netcdf_define_header( callmode, extend, av )
388
389#if defined( __netcdf )
390
391    USE arrays_3d,                                                             &
392        ONLY:  zu, zw
393
394    USE biometeorology_mod,                                                    &
395        ONLY:  bio_define_netcdf_grid
396
397    USE chemistry_model_mod,                                                   &
398        ONLY:  chem_define_netcdf_grid
399
400    USE basic_constants_and_equations_mod,                                     &
401        ONLY:  pi
402
403    USE control_parameters,                                                    &
404        ONLY:  agent_time_unlimited, air_chemistry, averaging_interval,        &
405               averaging_interval_pr, data_output_pr, domask, dopr_n,          &
406               dopr_time_count, dopts_time_count, dots_time_count,             &
407               do2d, do2d_at_begin, do2d_xz_time_count, do3d, do3d_at_begin,   &
408               do2d_yz_time_count, dt_data_output_av, dt_do2d_xy, dt_do2d_xz,  &
409               dt_do2d_yz, dt_do3d, dt_write_agent_data, mask_size,            &
410               do2d_xy_time_count, do3d_time_count, domask_time_count,         &
411               end_time, indoor_model, land_surface,                           &
412               mask_size_l, mask_i, mask_i_global, mask_j, mask_j_global,      &
413               mask_k_global, mask_surface,                                    &
414               message_string, ntdim_2d_xy, ntdim_2d_xz,                       &
415               ntdim_2d_yz, ntdim_3d, nz_do3d, ocean_mode, plant_canopy,       &
416               run_description_header, salsa, section, simulated_time,         &
417               simulated_time_at_begin, skip_time_data_output_av,              &
418               skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,        &
419               skip_time_do3d, topography, num_leg, num_var_fl,                &
420               urban_surface
421
422    USE diagnostic_output_quantities_mod,                                      &
423        ONLY:  doq_define_netcdf_grid
424
425    USE grid_variables,                                                        &
426        ONLY:  dx, dy, zu_s_inner, zw_w_inner
427
428    USE gust_mod,                                                              &
429        ONLY: gust_define_netcdf_grid, gust_module_enabled
430
431    USE indices,                                                               &
432        ONLY:  nx, nxl, nxr, ny, nys, nyn, nz ,nzb, nzt
433
434    USE kinds
435
436    USE indoor_model_mod,                                                      &
437        ONLY: im_define_netcdf_grid
438
439    USE land_surface_model_mod,                                                &
440        ONLY: lsm_define_netcdf_grid, nzb_soil, nzt_soil, nzs, zs
441
442    USE ocean_mod,                                                             &
443        ONLY:  ocean_define_netcdf_grid
444
445    USE pegrid
446
447    USE particle_attributes,                                                   &
448        ONLY:  number_of_particle_groups
449
450    USE plant_canopy_model_mod,                                                &
451        ONLY:  pch_index, pcm_define_netcdf_grid
452
453    USE profil_parameter,                                                      &
454        ONLY:  crmax, cross_profiles, dopr_index, profile_columns, profile_rows
455
456    USE radiation_model_mod,                                                   &
457        ONLY: radiation, radiation_define_netcdf_grid
458
459    USE salsa_mod,                                                             &
460        ONLY:  salsa_define_netcdf_grid
461
462    USE spectra_mod,                                                           &
463        ONLY:  averaging_interval_sp, comp_spectra_level, data_output_sp, dosp_time_count, spectra_direction
464
465    USE statistics,                                                            &
466        ONLY:  hom, statistic_regions
467
468    USE turbulence_closure_mod,                                                &
469        ONLY:  tcm_define_netcdf_grid
470
471    USE urban_surface_mod,                                                     &
472        ONLY:  usm_define_netcdf_grid
473
474    USE user,                                                                  &
475        ONLY:  user_module_enabled, user_define_netcdf_grid
476
477
478
479    IMPLICIT NONE
480
481    CHARACTER (LEN=3)              ::  suffix                !<
482    CHARACTER (LEN=2), INTENT (IN) ::  callmode              !<
483    CHARACTER (LEN=4)              ::  grid_x                !<
484    CHARACTER (LEN=4)              ::  grid_y                !<
485    CHARACTER (LEN=4)              ::  grid_z                !<
486    CHARACTER (LEN=6)              ::  mode                  !<
487    CHARACTER (LEN=10)             ::  precision             !<
488    CHARACTER (LEN=10)             ::  var                   !<
489    CHARACTER (LEN=20)             ::  netcdf_var_name       !<
490    CHARACTER (LEN=varnamelength)  ::  trimvar               !< TRIM of output-variable string
491    CHARACTER (LEN=80)             ::  time_average_text     !<
492    CHARACTER (LEN=4000)           ::  char_cross_profiles   !<
493    CHARACTER (LEN=4000)           ::  var_list              !<
494    CHARACTER (LEN=4000)           ::  var_list_old          !<
495
496    CHARACTER (LEN=100), DIMENSION(1:crmax) ::  cross_profiles_adj   !<
497    CHARACTER (LEN=100), DIMENSION(1:crmax) ::  cross_profiles_char  !<
498
499    INTEGER(iwp) ::  av                                      !<
500    INTEGER(iwp) ::  cross_profiles_count                    !<
501    INTEGER(iwp) ::  cross_profiles_maxi                     !<
502    INTEGER(iwp) ::  delim                                   !<
503    INTEGER(iwp) ::  delim_old                               !<
504    INTEGER(iwp) ::  file_id                                 !<
505    INTEGER(iwp) ::  i                                       !<
506    INTEGER(iwp) ::  id_last                                 !<
507    INTEGER(iwp) ::  id_x                                    !<
508    INTEGER(iwp) ::  id_y                                    !<
509    INTEGER(iwp) ::  id_z                                    !<
510    INTEGER(iwp) ::  j                                       !<
511    INTEGER(iwp) ::  k                                       !<
512    INTEGER(iwp) ::  kk                                      !<
513    INTEGER(iwp) ::  mid                                     !< masked output running index
514    INTEGER(iwp) ::  ns                                      !<
515    INTEGER(iwp) ::  ns_do                                   !< actual value of ns for soil model data
516    INTEGER(iwp) ::  ns_old                                  !<
517    INTEGER(iwp) ::  ntime_count                             !< number of time levels found in file
518    INTEGER(iwp) ::  nz_old                                  !<
519    INTEGER(iwp) ::  l                                       !<
520
521    INTEGER(iwp), SAVE ::  oldmode                           !<
522
523    INTEGER(iwp), DIMENSION(1) ::  id_dim_time_old           !<
524    INTEGER(iwp), DIMENSION(1) ::  id_dim_x_yz_old           !<
525    INTEGER(iwp), DIMENSION(1) ::  id_dim_y_xz_old           !<
526    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_sp_old          !<
527    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_xy_old          !<
528    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_3d_old          !<
529    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_mask_old        !<
530
531
532    INTEGER(iwp), DIMENSION(1:crmax) ::  cross_profiles_numb !<
533
534    LOGICAL ::  found                                        !<
535
536    LOGICAL, INTENT (INOUT) ::  extend                       !<
537
538    LOGICAL, SAVE ::  init_netcdf = .FALSE.                  !<
539
540    REAL(wp) ::  cos_ra                                      !< cosine of rotation_angle
541    REAL(wp) ::  eutm                                        !< easting (UTM)
542    REAL(wp) ::  nutm                                        !< northing (UTM)
543    REAL(wp) ::  shift_x                                     !< shift of x coordinate
544    REAL(wp) ::  shift_y                                     !< shift of y coordinate
545    REAL(wp) ::  sin_ra                                      !< sine of rotation_angle
546
547    REAL(wp), DIMENSION(1) ::  last_time_coordinate          !< last time value in file
548    REAL(wp), DIMENSION(8) ::  crs_list                      !< list of coord_ref_sys values
549
550    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  netcdf_data    !<
551    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  netcdf_data_2d !<
552    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lat            !< latitude
553    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lon            !< longitude
554
555
556!
557!-- Initializing actions
558    IF ( .NOT. init_netcdf )  THEN
559!
560!--    Check and set accuracy for netCDF output. First set default value
561       nc_precision = NF90_REAL4
562
563       i = 1
564       DO  WHILE ( netcdf_precision(i) /= ' ' )
565          j = INDEX( netcdf_precision(i), '_' )
566          IF ( j == 0 )  THEN
567             WRITE ( message_string, * ) 'netcdf_precision must contain a ', &
568                                         '"_"netcdf_precision(', i, ')="',   &
569                                         TRIM( netcdf_precision(i) ),'"'
570             CALL message( 'netcdf_define_header', 'PA0241', 2, 2, 0, 6, 0 )
571          ENDIF
572
573          var       = netcdf_precision(i)(1:j-1)
574          precision = netcdf_precision(i)(j+1:)
575
576          IF ( precision == 'NF90_REAL4' )  THEN
577             j = NF90_REAL4
578          ELSEIF ( precision == 'NF90_REAL8' )  THEN
579             j = NF90_REAL8
580          ELSE
581             WRITE ( message_string, * ) 'illegal netcdf precision: ',  &
582                                         'netcdf_precision(', i, ')="', &
583                                         TRIM( netcdf_precision(i) ),'"'
584             CALL message( 'netcdf_define_header', 'PA0242', 1, 2, 0, 6, 0 )
585          ENDIF
586
587          SELECT CASE ( var )
588             CASE ( 'xy' )
589                nc_precision(1) = j
590             CASE ( 'xz' )
591                nc_precision(2) = j
592             CASE ( 'yz' )
593                nc_precision(3) = j
594             CASE ( '2d' )
595                nc_precision(1:3) = j
596             CASE ( '3d' )
597                nc_precision(4) = j
598             CASE ( 'pr' )
599                nc_precision(5) = j
600             CASE ( 'ts' )
601                nc_precision(6) = j
602             CASE ( 'sp' )
603                nc_precision(7) = j
604             CASE ( 'prt' )
605                nc_precision(8) = j
606             CASE ( 'masks' )
607                nc_precision(11) = j
608             CASE ( 'fl' )
609                nc_precision(9) = j
610             CASE ( 'all' )
611                nc_precision    = j
612
613             CASE DEFAULT
614                WRITE ( message_string, * ) 'unknown variable in ' //          &
615                                  'initialization_parameters ',                &
616                                  'assignment: netcdf_precision(', i, ')="',   &
617                                            TRIM( netcdf_precision(i) ),'"'
618                CALL message( 'netcdf_define_header', 'PA0243', 1, 2, 0, 6, 0 )
619
620          END SELECT
621
622          i = i + 1
623          IF ( i > 50 )  EXIT
624       ENDDO
625
626!
627!--    Check for allowed parameter range
628       IF ( netcdf_deflate < 0  .OR.  netcdf_deflate > 9 )  THEN
629          WRITE ( message_string, '(A,I3,A)' ) 'netcdf_deflate out of ' //     &
630                                      'range & given value: ', netcdf_deflate, &
631                                      ', allowed range: 0-9'
632          CALL message( 'netcdf_define_header', 'PA0355', 2, 2, 0, 6, 0 )
633       ENDIF
634!
635!--    Data compression does not work with parallel NetCDF/HDF5
636       IF ( netcdf_deflate > 0  .AND.  netcdf_data_format /= 3 )  THEN
637          message_string = 'netcdf_deflate reset to 0'
638          CALL message( 'netcdf_define_header', 'PA0356', 0, 1, 0, 6, 0 )
639
640          netcdf_deflate = 0
641       ENDIF
642
643       init_netcdf = .TRUE.
644
645    ENDIF
646!
647!-- Convert coord_ref_sys into vector (used for lat/lon calculation)
648    crs_list = (/ coord_ref_sys%semi_major_axis,                  &
649                  coord_ref_sys%inverse_flattening,               &
650                  coord_ref_sys%longitude_of_prime_meridian,      &
651                  coord_ref_sys%longitude_of_central_meridian,    &
652                  coord_ref_sys%scale_factor_at_central_meridian, &
653                  coord_ref_sys%latitude_of_projection_origin,    &
654                  coord_ref_sys%false_easting,                    &
655                  coord_ref_sys%false_northing /)
656
657!
658!-- Determine the mode to be processed
659    IF ( extend )  THEN
660       mode = callmode // '_ext'
661    ELSE
662       mode = callmode // '_new'
663    ENDIF
664
665!
666!-- Select the mode to be processed. Possibilities are 3d, ma (mask), xy, xz,
667!-- yz, pr (profiles), ps (particle timeseries), fl (flight data), ts
668!-- (timeseries) or sp (spectra)
669    SELECT CASE ( mode )
670
671       CASE ( 'ma_new' )
672
673!
674!--       decompose actual parameter file_id (=formal parameter av) into
675!--       mid and av
676          file_id = av
677          IF ( file_id <= 200+max_masks )  THEN
678             mid = file_id - 200
679             av = 0
680          ELSE
681             mid = file_id - (200+max_masks)
682             av = 1
683          ENDIF
684
685!
686!--       Define some global attributes of the dataset
687          IF ( av == 0 )  THEN
688             CALL netcdf_create_global_atts( id_set_mask(mid,av), 'podsmasked', TRIM( run_description_header ), 464 )
689             time_average_text = ' '
690          ELSE
691             CALL netcdf_create_global_atts( id_set_mask(mid,av), 'podsmasked', TRIM( run_description_header ), 464 )
692             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
693             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'time_avg',   &
694                                     TRIM( time_average_text ) )
695             CALL netcdf_handle_error( 'netcdf_define_header', 466 )
696          ENDIF
697
698!
699!--       Define time coordinate for volume data (unlimited dimension)
700          CALL netcdf_create_dim( id_set_mask(mid,av), 'time', NF90_UNLIMITED, &
701                                  id_dim_time_mask(mid,av), 467 )
702          CALL netcdf_create_var( id_set_mask(mid,av),                         &
703                                  (/ id_dim_time_mask(mid,av) /), 'time',      &
704                                  NF90_DOUBLE, id_var_time_mask(mid,av),       &
705                                 'seconds', 'time', 468, 469, 000 )
706          CALL netcdf_create_att( id_set_mask(mid,av), id_var_time_mask(mid,av), 'standard_name', 'time', 000)
707          CALL netcdf_create_att( id_set_mask(mid,av), id_var_time_mask(mid,av), 'axis', 'T', 000)
708
709!
710!--       Define spatial dimensions and coordinates:
711          IF ( mask_surface(mid) )  THEN
712!
713!--          In case of terrain-following output, the vertical dimensions are
714!--          indices, not meters
715             CALL netcdf_create_dim( id_set_mask(mid,av), 'ku_above_surf',     &
716                                     mask_size(mid,3), id_dim_zu_mask(mid,av), &
717                                     470 )
718             CALL netcdf_create_var( id_set_mask(mid,av),                      &
719                                     (/ id_dim_zu_mask(mid,av) /),             &
720                                     'ku_above_surf',                          &
721                                     NF90_DOUBLE, id_var_zu_mask(mid,av),      &
722                                     '1', 'grid point above terrain',          &
723                                     471, 472, 000 )
724             CALL netcdf_create_dim( id_set_mask(mid,av), 'kw_above_surf',     &
725                                     mask_size(mid,3), id_dim_zw_mask(mid,av), &
726                                     473 )
727             CALL netcdf_create_var( id_set_mask(mid,av),                      &
728                                     (/ id_dim_zw_mask(mid,av) /),             &
729                                     'kw_above_surf',                          &
730                                     NF90_DOUBLE, id_var_zw_mask(mid,av),      &
731                                    '1', 'grid point above terrain',           &
732                                    474, 475, 000 )
733          ELSE
734!
735!--          Define vertical coordinate grid (zu grid)
736             CALL netcdf_create_dim( id_set_mask(mid,av), 'zu_3d',             &
737                                     mask_size(mid,3), id_dim_zu_mask(mid,av), &
738                                     470 )
739             CALL netcdf_create_var( id_set_mask(mid,av),                      &
740                                     (/ id_dim_zu_mask(mid,av) /), 'zu_3d',    &
741                                     NF90_DOUBLE, id_var_zu_mask(mid,av),      &
742                                     'meters', '', 471, 472, 000 )
743!
744!--          Define vertical coordinate grid (zw grid)
745             CALL netcdf_create_dim( id_set_mask(mid,av), 'zw_3d',             &
746                                     mask_size(mid,3), id_dim_zw_mask(mid,av), &
747                                     473 )
748             CALL netcdf_create_var( id_set_mask(mid,av),                      &
749                                     (/ id_dim_zw_mask(mid,av) /), 'zw_3d',    &
750                                     NF90_DOUBLE, id_var_zw_mask(mid,av),      &
751                                    'meters', '', 474, 475, 000 )
752          ENDIF
753!
754!--       Define x-axis (for scalar position)
755          CALL netcdf_create_dim( id_set_mask(mid,av), 'x', mask_size(mid,1),  &
756                                  id_dim_x_mask(mid,av), 476 )
757          CALL netcdf_create_var( id_set_mask(mid,av),                         &
758                                  (/ id_dim_x_mask(mid,av) /), 'x',            &
759                                  NF90_DOUBLE, id_var_x_mask(mid,av),          &
760                                  'meters', '', 477, 478, 000 )
761!
762!--       Define x-axis (for u position)
763          CALL netcdf_create_dim( id_set_mask(mid,av), 'xu', mask_size(mid,1), &
764                                  id_dim_xu_mask(mid,av), 479 )
765          CALL netcdf_create_var( id_set_mask(mid,av),                         &
766                                  (/ id_dim_xu_mask(mid,av) /), 'xu',          &
767                                  NF90_DOUBLE, id_var_xu_mask(mid,av),         &
768                                  'meters', '', 480, 481, 000 )
769!
770!--       Define y-axis (for scalar position)
771          CALL netcdf_create_dim( id_set_mask(mid,av), 'y', mask_size(mid,2),  &
772                                  id_dim_y_mask(mid,av), 482 )
773          CALL netcdf_create_var( id_set_mask(mid,av),                         &
774                                  (/ id_dim_y_mask(mid,av) /), 'y',            &
775                                  NF90_DOUBLE, id_var_y_mask(mid,av),          &
776                                  'meters', '', 483, 484, 000 )
777!
778!--       Define y-axis (for v position)
779          CALL netcdf_create_dim( id_set_mask(mid,av), 'yv', mask_size(mid,2), &
780                                  id_dim_yv_mask(mid,av), 485 )
781          CALL netcdf_create_var( id_set_mask(mid,av),                         &
782                                  (/ id_dim_yv_mask(mid,av) /),                &
783                                  'yv', NF90_DOUBLE, id_var_yv_mask(mid,av),   &
784                                  'meters', '', 486, 487, 000 )
785!
786!--       Define UTM and geographic coordinates
787          CALL define_geo_coordinates( id_set_mask(mid,av),               &
788                  (/ id_dim_x_mask(mid,av), id_dim_xu_mask(mid,av) /),    &
789                  (/ id_dim_y_mask(mid,av), id_dim_yv_mask(mid,av) /),    &
790                  id_var_eutm_mask(:,mid,av), id_var_nutm_mask(:,mid,av), &
791                  id_var_lat_mask(:,mid,av), id_var_lon_mask(:,mid,av)    )
792!
793!--       Define coordinate-reference system
794          CALL netcdf_create_crs( id_set_mask(mid,av), 000 )
795!
796!--       In case of non-flat topography define 2d-arrays containing the height
797!--       information. Only for parallel netcdf output.
798          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
799               netcdf_data_format > 4 )  THEN
800!
801!--          Define zusi = zu(nzb_s_inner)
802             CALL netcdf_create_var( id_set_mask(mid,av),                      &
803                                     (/ id_dim_x_mask(mid,av),                 &
804                                        id_dim_y_mask(mid,av) /), 'zusi',      &
805                                     NF90_DOUBLE, id_var_zusi_mask(mid,av),    &
806                                     'meters', 'zu(nzb_s_inner)', 488, 489,    &
807                                     490 )
808!
809!--          Define zwwi = zw(nzb_w_inner)
810             CALL netcdf_create_var( id_set_mask(mid,av),                      &
811                                     (/ id_dim_x_mask(mid,av),                 &
812                                        id_dim_y_mask(mid,av) /), 'zwwi',      &
813                                     NF90_DOUBLE, id_var_zwwi_mask(mid,av),    &
814                                     'meters', 'zw(nzb_w_inner)', 491, 492,    &
815                                     493 )
816          ENDIF
817
818          IF ( land_surface )  THEN
819!
820!--          Define vertical coordinate grid (zw grid)
821             CALL netcdf_create_dim( id_set_mask(mid,av), 'zs_3d',             &
822                                     mask_size(mid,3), id_dim_zs_mask(mid,av), &
823                                     536 )
824             CALL netcdf_create_var( id_set_mask(mid,av),                      &
825                                     (/ id_dim_zs_mask(mid,av) /), 'zs_3d',    &
826                                     NF90_DOUBLE, id_var_zs_mask(mid,av),      &
827                                     'meters', '', 537, 555, 000 )
828          ENDIF
829
830!
831!--       Define the variables
832          var_list = ';'
833          i = 1
834
835          DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
836
837             trimvar = TRIM( domask(mid,av,i) )
838!
839!--          Check for the grid
840             found = .FALSE.
841             SELECT CASE ( trimvar )
842!
843!--             Most variables are defined on the scalar grid
844                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',                &
845                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
846                       's', 'theta', 'thetal', 'thetav' )
847
848                   grid_x = 'x'
849                   grid_y = 'y'
850                   grid_z = 'zu'
851!
852!--             u grid
853                CASE ( 'u' )
854
855                   grid_x = 'xu'
856                   grid_y = 'y'
857                   grid_z = 'zu'
858!
859!--             v grid
860                CASE ( 'v' )
861
862                   grid_x = 'x'
863                   grid_y = 'yv'
864                   grid_z = 'zu'
865!
866!--             w grid
867                CASE ( 'w' )
868
869                   grid_x = 'x'
870                   grid_y = 'y'
871                   grid_z = 'zw'
872
873
874                CASE DEFAULT
875!
876!--                Check for quantities defined in other modules
877                   CALL tcm_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
878
879                   IF ( .NOT. found  .AND.  air_chemistry )  THEN
880                      CALL chem_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
881                   ENDIF
882
883                   IF ( .NOT. found )                                          &
884                      CALL doq_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
885
886                   IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
887                      CALL gust_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
888                   ENDIF
889
890                   IF ( .NOT. found  .AND. land_surface )  THEN
891                      CALL lsm_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
892                   ENDIF
893
894                   IF ( .NOT. found  .AND.  ocean_mode )  THEN
895                      CALL ocean_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
896                   ENDIF
897
898                   IF ( .NOT. found  .AND.  plant_canopy )  THEN
899                      CALL pcm_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
900                   ENDIF
901
902                   IF ( .NOT. found  .AND.  radiation )  THEN
903                      CALL radiation_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
904                   ENDIF
905!
906!--                Check for SALSA quantities
907                   IF ( .NOT. found  .AND.  salsa )  THEN
908                      CALL salsa_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
909                   ENDIF
910
911                   IF ( .NOT. found  .AND.  urban_surface )  THEN
912                      CALL usm_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
913                   ENDIF
914!
915!--                Now check for user-defined quantities
916                   IF ( .NOT. found  .AND.  user_module_enabled )  THEN
917                      CALL user_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
918                   ENDIF
919
920                   IF ( .NOT. found )  THEN
921                      WRITE ( message_string, * ) 'no grid defined for variable ', TRIM( trimvar )
922                      CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, 6, 0 )
923                   ENDIF
924
925             END SELECT
926
927!
928!--          Select the respective dimension ids
929             IF ( grid_x == 'x' )  THEN
930                id_x = id_dim_x_mask(mid,av)
931             ELSEIF ( grid_x == 'xu' )  THEN
932                id_x = id_dim_xu_mask(mid,av)
933             ENDIF
934
935             IF ( grid_y == 'y' )  THEN
936                id_y = id_dim_y_mask(mid,av)
937             ELSEIF ( grid_y == 'yv' )  THEN
938                id_y = id_dim_yv_mask(mid,av)
939             ENDIF
940
941             IF ( grid_z == 'zu' )  THEN
942                id_z = id_dim_zu_mask(mid,av)
943             ELSEIF ( grid_z == 'zw' )  THEN
944                id_z = id_dim_zw_mask(mid,av)
945             ELSEIF ( grid_z == "zs" )  THEN
946                id_z = id_dim_zs_mask(mid,av)
947             ENDIF
948
949!
950!--          Define the grid
951             CALL netcdf_create_var( id_set_mask(mid,av), (/ id_x, id_y, id_z, &
952                                     id_dim_time_mask(mid,av) /),              &
953                                     domask(mid,av,i), nc_precision(11),       &
954                                     id_var_domask(mid,av,i),                  &
955                                     TRIM( domask_unit(mid,av,i) ),            &
956                                     domask(mid,av,i), 494, 495, 496, .TRUE. )
957
958             var_list = TRIM( var_list ) // TRIM( domask(mid,av,i) ) // ';'
959
960             i = i + 1
961
962          ENDDO
963
964!
965!--       No arrays to output
966          IF ( i == 1 )  RETURN
967
968!
969!--       Write the list of variables as global attribute (this is used by
970!--       restart runs and by combine_plot_fields)
971          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
972                                  'VAR_LIST', var_list )
973          CALL netcdf_handle_error( 'netcdf_define_header', 497 )
974
975!
976!--       Leave netCDF define mode
977          nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
978          CALL netcdf_handle_error( 'netcdf_define_header', 498 )
979
980!
981!--       Write data for x (shifted by +dx/2) and xu axis
982          ALLOCATE( netcdf_data(mask_size(mid,1)) )
983
984          netcdf_data = ( mask_i_global(mid,:mask_size(mid,1)) + 0.5_wp ) * dx
985
986          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_x_mask(mid,av), &
987                                  netcdf_data, start = (/ 1 /),               &
988                                  count = (/ mask_size(mid,1) /) )
989          CALL netcdf_handle_error( 'netcdf_define_header', 499 )
990
991          netcdf_data = mask_i_global(mid,:mask_size(mid,1)) * dx
992
993          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_xu_mask(mid,av),&
994                                  netcdf_data, start = (/ 1 /),               &
995                                  count = (/ mask_size(mid,1) /) )
996          CALL netcdf_handle_error( 'netcdf_define_header', 500 )
997
998          DEALLOCATE( netcdf_data )
999
1000!
1001!--       Write data for y (shifted by +dy/2) and yv axis
1002          ALLOCATE( netcdf_data(mask_size(mid,2)) )
1003
1004          netcdf_data = ( mask_j_global(mid,:mask_size(mid,2)) + 0.5_wp ) * dy
1005
1006          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_y_mask(mid,av), &
1007                                  netcdf_data, start = (/ 1 /),               &
1008                                  count = (/ mask_size(mid,2) /))
1009          CALL netcdf_handle_error( 'netcdf_define_header', 501 )
1010
1011          netcdf_data = mask_j_global(mid,:mask_size(mid,2)) * dy
1012
1013          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_yv_mask(mid,av), &
1014                                  netcdf_data, start = (/ 1 /),    &
1015                                  count = (/ mask_size(mid,2) /))
1016          CALL netcdf_handle_error( 'netcdf_define_header', 502 )
1017
1018          DEALLOCATE( netcdf_data )
1019
1020!
1021!--       Write UTM coordinates
1022          IF ( init_model%rotation_angle == 0.0_wp )  THEN
1023!
1024!--          1D in case of no rotation
1025             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
1026!
1027!--          x coordinates
1028             ALLOCATE( netcdf_data(mask_size(mid,1)) )
1029             DO  k = 0, 2
1030!
1031!--             Scalar grid points
1032                IF ( k == 0 )  THEN
1033                   shift_x = 0.5
1034!
1035!--             u grid points
1036                ELSEIF ( k == 1 )  THEN
1037                   shift_x = 0.0
1038!
1039!--             v grid points
1040                ELSEIF ( k == 2 )  THEN
1041                   shift_x = 0.5
1042                ENDIF
1043
1044                netcdf_data = init_model%origin_x + cos_ra                     &
1045                       * ( mask_i_global(mid,:mask_size(mid,1)) + shift_x ) * dx
1046
1047                nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
1048                                        id_var_eutm_mask(k,mid,av), &
1049                                        netcdf_data, start = (/ 1 /), &
1050                                        count = (/ mask_size(mid,1) /) )
1051                CALL netcdf_handle_error( 'netcdf_define_header', 555 )
1052
1053             ENDDO
1054             DEALLOCATE( netcdf_data )
1055!
1056!--          y coordinates
1057             ALLOCATE( netcdf_data(mask_size(mid,2)) )
1058             DO  k = 0, 2
1059!
1060!--             Scalar grid points
1061                IF ( k == 0 )  THEN
1062                   shift_y = 0.5
1063!
1064!--             u grid points
1065                ELSEIF ( k == 1 )  THEN
1066                   shift_y = 0.5
1067!
1068!--             v grid points
1069                ELSEIF ( k == 2 )  THEN
1070                   shift_y = 0.0
1071                ENDIF
1072
1073                netcdf_data = init_model%origin_y + cos_ra                     &
1074                       * ( mask_j_global(mid,:mask_size(mid,2)) + shift_y ) * dy
1075
1076                nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
1077                                        id_var_nutm_mask(k,mid,av), &
1078                                        netcdf_data, start = (/ 1 /), &
1079                                        count = (/ mask_size(mid,2) /) )
1080                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
1081
1082             ENDDO
1083             DEALLOCATE( netcdf_data )
1084
1085          ELSE
1086!
1087!--          2D in case of rotation
1088             ALLOCATE( netcdf_data_2d(mask_size(mid,1),mask_size(mid,2)) )
1089             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
1090             sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
1091
1092             DO  k = 0, 2
1093!
1094!--            Scalar grid points
1095               IF ( k == 0 )  THEN
1096                  shift_x = 0.5 ; shift_y = 0.5
1097!
1098!--            u grid points
1099               ELSEIF ( k == 1 )  THEN
1100                  shift_x = 0.0 ; shift_y = 0.5
1101!
1102!--            v grid points
1103               ELSEIF ( k == 2 )  THEN
1104                  shift_x = 0.5 ; shift_y = 0.0
1105               ENDIF
1106
1107               DO  j = 1, mask_size(mid,2)
1108                  DO  i = 1, mask_size(mid,1)
1109                     netcdf_data_2d(i,j) = init_model%origin_x                 &
1110                           + cos_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
1111                           + sin_ra * ( mask_j_global(mid,j) + shift_y ) * dy
1112                  ENDDO
1113               ENDDO
1114
1115               nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
1116                                       id_var_eutm_mask(k,mid,av), &
1117                                       netcdf_data_2d, start = (/ 1, 1 /), &
1118                                       count = (/ mask_size(mid,1), &
1119                                                  mask_size(mid,2) /) )
1120               CALL netcdf_handle_error( 'netcdf_define_header', 555 )
1121
1122               DO  j = 1, mask_size(mid,2)
1123                  DO  i = 1, mask_size(mid,1)
1124                     netcdf_data_2d(i,j) = init_model%origin_y                 &
1125                           - sin_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
1126                           + cos_ra * ( mask_j_global(mid,j) + shift_y ) * dy
1127                  ENDDO
1128               ENDDO
1129
1130               nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
1131                                       id_var_nutm_mask(k,mid,av), &
1132                                       netcdf_data_2d, start = (/ 1, 1 /), &
1133                                       count = (/ mask_size(mid,1), &
1134                                                  mask_size(mid,2) /) )
1135               CALL netcdf_handle_error( 'netcdf_define_header', 556 )
1136
1137             ENDDO
1138             DEALLOCATE( netcdf_data_2d )
1139          ENDIF
1140!
1141!--       Write lon and lat data.
1142          ALLOCATE( lat(mask_size(mid,1),mask_size(mid,2)) )
1143          ALLOCATE( lon(mask_size(mid,1),mask_size(mid,2)) )
1144          cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
1145          sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
1146
1147          DO  k = 0, 2
1148!
1149!--          Scalar grid points
1150             IF ( k == 0 )  THEN
1151                shift_x = 0.5 ; shift_y = 0.5
1152!
1153!--          u grid points
1154             ELSEIF ( k == 1 )  THEN
1155                shift_x = 0.0 ; shift_y = 0.5
1156!
1157!--          v grid points
1158             ELSEIF ( k == 2 )  THEN
1159                shift_x = 0.5 ; shift_y = 0.0
1160             ENDIF
1161
1162             DO  j = 1, mask_size(mid,2)
1163                DO  i = 1, mask_size(mid,1)
1164                   eutm = init_model%origin_x                               &
1165                        + cos_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
1166                        + sin_ra * ( mask_j_global(mid,j) + shift_y ) * dy
1167                   nutm = init_model%origin_y                               &
1168                        - sin_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
1169                        + cos_ra * ( mask_j_global(mid,j) + shift_y ) * dy
1170
1171                   CALL  convert_utm_to_geographic( crs_list,          &
1172                                                    eutm, nutm,        &
1173                                                    lon(i,j), lat(i,j) )
1174                ENDDO
1175             ENDDO
1176
1177             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),           &
1178                                     id_var_lon_mask(k,mid,av),     &
1179                                     lon, start = (/ 1, 1 /),       &
1180                                     count = (/ mask_size(mid,1),   &
1181                                                mask_size(mid,2) /) )
1182             CALL netcdf_handle_error( 'netcdf_define_header', 556 )
1183
1184             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),           &
1185                                     id_var_lat_mask(k,mid,av),     &
1186                                     lat, start = (/ 1, 1 /),       &
1187                                     count = (/ mask_size(mid,1),   &
1188                                                mask_size(mid,2) /) )
1189             CALL netcdf_handle_error( 'netcdf_define_header', 556 )
1190          ENDDO
1191
1192          DEALLOCATE( lat )
1193          DEALLOCATE( lon )
1194!
1195!--       Write zu and zw data (vertical axes)
1196          ALLOCATE( netcdf_data(mask_size(mid,3)) )
1197
1198          IF ( mask_surface(mid) )  THEN
1199
1200             netcdf_data = mask_k_global(mid,:mask_size(mid,3))
1201
1202             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
1203                                     netcdf_data, start = (/ 1 /), &
1204                                     count = (/ mask_size(mid,3) /) )
1205             CALL netcdf_handle_error( 'netcdf_define_header', 503 )
1206
1207             netcdf_data = mask_k_global(mid,:mask_size(mid,3))
1208
1209             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
1210                                     netcdf_data, start = (/ 1 /), &
1211                                     count = (/ mask_size(mid,3) /) )
1212             CALL netcdf_handle_error( 'netcdf_define_header', 504 )
1213
1214          ELSE
1215
1216             netcdf_data = zu( mask_k_global(mid,:mask_size(mid,3)) )
1217
1218             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
1219                                     netcdf_data, start = (/ 1 /), &
1220                                     count = (/ mask_size(mid,3) /) )
1221             CALL netcdf_handle_error( 'netcdf_define_header', 503 )
1222
1223             netcdf_data = zw( mask_k_global(mid,:mask_size(mid,3)) )
1224
1225             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
1226                                     netcdf_data, start = (/ 1 /), &
1227                                     count = (/ mask_size(mid,3) /) )
1228             CALL netcdf_handle_error( 'netcdf_define_header', 504 )
1229
1230          ENDIF
1231
1232          DEALLOCATE( netcdf_data )
1233
1234!
1235!--       In case of non-flat topography write height information
1236          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
1237               netcdf_data_format > 4 )  THEN
1238
1239             ALLOCATE( netcdf_data_2d(mask_size_l(mid,1),mask_size_l(mid,2)) )
1240             netcdf_data_2d = zu_s_inner( mask_i(mid,:mask_size_l(mid,1)),     &
1241                                          mask_j(mid,:mask_size_l(mid,2)) )
1242
1243             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),                      &
1244                                     id_var_zusi_mask(mid,av),                 &
1245                                     netcdf_data_2d,                           &
1246                                     start = (/ 1, 1 /),                       &
1247                                     count = (/ mask_size_l(mid,1),            &
1248                                                mask_size_l(mid,2) /) )
1249             CALL netcdf_handle_error( 'netcdf_define_header', 505 )
1250
1251             netcdf_data_2d = zw_w_inner( mask_i(mid,:mask_size_l(mid,1)),     &
1252                                          mask_j(mid,:mask_size_l(mid,2)) )
1253
1254             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),                      &
1255                                     id_var_zwwi_mask(mid,av),                 &
1256                                     netcdf_data_2d,                           &
1257                                     start = (/ 1, 1 /),                       &
1258                                     count = (/ mask_size_l(mid,1),            &
1259                                                mask_size_l(mid,2) /) )
1260             CALL netcdf_handle_error( 'netcdf_define_header', 506 )
1261
1262             DEALLOCATE( netcdf_data_2d )
1263
1264          ENDIF
1265!
1266!--       soil is not in masked output for now - disable temporary this block
1267!          IF ( land_surface )  THEN
1268!
1269!--          Write zs data (vertical axes for soil model), use negative values
1270!--          to indicate soil depth
1271!             ALLOCATE( netcdf_data(mask_size(mid,3)) )
1272!
1273!             netcdf_data = zs( mask_k_global(mid,:mask_size(mid,3)) )
1274!
1275!             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zs_mask(mid,av), &
1276!                                     netcdf_data, start = (/ 1 /), &
1277!                                     count = (/ mask_size(mid,3) /) )
1278!             CALL netcdf_handle_error( 'netcdf_define_header', 538 )
1279!
1280!             DEALLOCATE( netcdf_data )
1281!
1282!          ENDIF
1283
1284!
1285!--       restore original parameter file_id (=formal parameter av) into av
1286          av = file_id
1287
1288
1289       CASE ( 'ma_ext' )
1290
1291!
1292!--       decompose actual parameter file_id (=formal parameter av) into
1293!--       mid and av
1294          file_id = av
1295          IF ( file_id <= 200+max_masks )  THEN
1296             mid = file_id - 200
1297             av = 0
1298          ELSE
1299             mid = file_id - (200+max_masks)
1300             av = 1
1301          ENDIF
1302
1303!
1304!--       Get the list of variables and compare with the actual run.
1305!--       First var_list_old has to be reset, since GET_ATT does not assign
1306!--       trailing blanks.
1307          var_list_old = ' '
1308          nc_stat = NF90_GET_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'VAR_LIST',&
1309                                  var_list_old )
1310          CALL netcdf_handle_error( 'netcdf_define_header', 507 )
1311
1312          var_list = ';'
1313          i = 1
1314          DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
1315             var_list = TRIM(var_list) // TRIM( domask(mid,av,i) ) // ';'
1316             i = i + 1
1317          ENDDO
1318
1319          IF ( av == 0 )  THEN
1320             var = '(mask)'
1321          ELSE
1322             var = '(mask_av)'
1323          ENDIF
1324
1325          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
1326             WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ),       &
1327                  ' data for mask', mid, ' from previous run found,',           &
1328                  '&but this file cannot be extended due to variable ',         &
1329                  'mismatch.&New file is created instead.'
1330             CALL message( 'define_netcdf_header', 'PA0335', 0, 1, 0, 6, 0 )
1331             extend = .FALSE.
1332             RETURN
1333          ENDIF
1334
1335!
1336!--       Get and compare the number of vertical gridpoints
1337          nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'zu_3d', &
1338                                    id_var_zu_mask(mid,av) )
1339          CALL netcdf_handle_error( 'netcdf_define_header', 508 )
1340
1341          nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av),     &
1342                                           id_var_zu_mask(mid,av),  &
1343                                           dimids = id_dim_zu_mask_old )
1344          CALL netcdf_handle_error( 'netcdf_define_header', 509 )
1345          id_dim_zu_mask(mid,av) = id_dim_zu_mask_old(1)
1346
1347          nc_stat = NF90_INQUIRE_DIMENSION( id_set_mask(mid,av),               &
1348                                            id_dim_zu_mask(mid,av),            &
1349                                            len = nz_old )
1350          CALL netcdf_handle_error( 'netcdf_define_header', 510 )
1351
1352          IF ( mask_size(mid,3) /= nz_old )  THEN
1353             WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ),      &
1354                  '&data for mask', mid, ' from previous run found,',          &
1355                  ' but this file cannot be extended due to mismatch in ',     &
1356                  ' number of vertical grid points.',                          &
1357                  '&New file is created instead.'
1358             CALL message( 'define_netcdf_header', 'PA0336', 0, 1, 0, 6, 0 )
1359             extend = .FALSE.
1360             RETURN
1361          ENDIF
1362
1363!
1364!--       Get the id of the time coordinate (unlimited coordinate) and its
1365!--       last index on the file. The next time level is plmask..count+1.
1366!--       The current time must be larger than the last output time
1367!--       on the file.
1368          nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'time',               &
1369                                    id_var_time_mask(mid,av) )
1370          CALL netcdf_handle_error( 'netcdf_define_header', 511 )
1371
1372          nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av),                &
1373                                           id_var_time_mask(mid,av),           &
1374                                           dimids = id_dim_time_old )
1375          CALL netcdf_handle_error( 'netcdf_define_header', 512 )
1376          id_dim_time_mask(mid,av) = id_dim_time_old(1)
1377
1378          nc_stat = NF90_INQUIRE_DIMENSION( id_set_mask(mid,av),               &
1379                                            id_dim_time_mask(mid,av),          &
1380                                            len = domask_time_count(mid,av) )
1381          CALL netcdf_handle_error( 'netcdf_define_header', 513 )
1382
1383          nc_stat = NF90_GET_VAR( id_set_mask(mid,av),                         &
1384                                  id_var_time_mask(mid,av),                    &
1385                                  last_time_coordinate,                        &
1386                                  start = (/ domask_time_count(mid,av) /),     &
1387                                  count = (/ 1 /) )
1388          CALL netcdf_handle_error( 'netcdf_define_header', 514 )
1389
1390          IF ( last_time_coordinate(1) >= simulated_time )  THEN
1391             WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ),      &
1392                  ' data for mask', mid, ' from previous run found,',          &
1393                  '&but this file cannot be extended because the current ',    &
1394                  'output time is less or equal than the last output time ',   &
1395                  'on this file.&New file is created instead.'
1396             CALL message( 'define_netcdf_header', 'PA0337', 0, 1, 0, 6, 0 )
1397             domask_time_count(mid,av) = 0
1398             extend = .FALSE.
1399             RETURN
1400          ENDIF
1401
1402!
1403!--       Dataset seems to be extendable.
1404!--       Now get the variable ids.
1405          i = 1
1406          DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
1407             nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), &
1408                                       TRIM( domask(mid,av,i) ), &
1409                                       id_var_domask(mid,av,i) )
1410             CALL netcdf_handle_error( 'netcdf_define_header', 515 )
1411             i = i + 1
1412          ENDDO
1413
1414!
1415!--       Update the title attribute on file
1416!--       In order to avoid 'data mode' errors if updated attributes are larger
1417!--       than their original size, NF90_PUT_ATT is called in 'define mode'
1418!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
1419!--       performance loss due to data copying; an alternative strategy would be
1420!--       to ensure equal attribute size in a job chain. Maybe revise later.
1421          IF ( av == 0 )  THEN
1422             time_average_text = ' '
1423          ELSE
1424             WRITE (time_average_text, '('', '',F7.1,'' s average'')')         &
1425                                                            averaging_interval
1426          ENDIF
1427          nc_stat = NF90_REDEF( id_set_mask(mid,av) )
1428          CALL netcdf_handle_error( 'netcdf_define_header', 516 )
1429          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'title',   &
1430                                  TRIM( run_description_header ) //            &
1431                                  TRIM( time_average_text ) )
1432          CALL netcdf_handle_error( 'netcdf_define_header', 517 )
1433          nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
1434          CALL netcdf_handle_error( 'netcdf_define_header', 518 )
1435          WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ),         &
1436               ' data for mask', mid, ' from previous run found.',             &
1437               ' &This file will be extended.'
1438          CALL message( 'define_netcdf_header', 'PA0338', 0, 0, 0, 6, 0 )
1439!
1440!--       restore original parameter file_id (=formal parameter av) into av
1441          av = file_id
1442
1443
1444       CASE ( '3d_new' )
1445
1446!
1447!--       Define some global attributes of the dataset
1448          IF ( av == 0 )  THEN
1449             CALL netcdf_create_global_atts( id_set_3d(av), '3d', TRIM( run_description_header ), 62 )
1450             time_average_text = ' '
1451          ELSE
1452             CALL netcdf_create_global_atts( id_set_3d(av), '3d_av', TRIM( run_description_header ), 62 )
1453             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
1454             nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'time_avg',   &
1455                                     TRIM( time_average_text ) )
1456             CALL netcdf_handle_error( 'netcdf_define_header', 63 )
1457          ENDIF
1458
1459!
1460!--       Define time coordinate for volume data.
1461!--       For parallel output the time dimensions has to be limited, otherwise
1462!--       the performance drops significantly.
1463          IF ( netcdf_data_format < 5 )  THEN
1464             CALL netcdf_create_dim( id_set_3d(av), 'time', NF90_UNLIMITED,    &
1465                                     id_dim_time_3d(av), 64 )
1466          ELSE
1467             CALL netcdf_create_dim( id_set_3d(av), 'time', ntdim_3d(av),      &
1468                                     id_dim_time_3d(av), 523 )
1469          ENDIF
1470
1471          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_time_3d(av) /),     &
1472                                  'time', NF90_DOUBLE, id_var_time_3d(av),     &
1473                                  'seconds', 'time', 65, 66, 00 )
1474          CALL netcdf_create_att( id_set_3d(av), id_var_time_3d(av), 'standard_name', 'time', 000)
1475          CALL netcdf_create_att( id_set_3d(av), id_var_time_3d(av), 'axis', 'T', 000)
1476!
1477!--       Define spatial dimensions and coordinates:
1478!--       Define vertical coordinate grid (zu grid)
1479          CALL netcdf_create_dim( id_set_3d(av), 'zu_3d', nz_do3d-nzb+1,       &
1480                                  id_dim_zu_3d(av), 67 )
1481          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zu_3d(av) /),       &
1482                                  'zu_3d', NF90_DOUBLE, id_var_zu_3d(av),      &
1483                                  'meters', '', 68, 69, 00 )
1484!
1485!--       Define vertical coordinate grid (zw grid)
1486          CALL netcdf_create_dim( id_set_3d(av), 'zw_3d', nz_do3d-nzb+1,       &
1487                                  id_dim_zw_3d(av), 70 )
1488          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zw_3d(av) /),       &
1489                                  'zw_3d', NF90_DOUBLE, id_var_zw_3d(av),      &
1490                                  'meters', '', 71, 72, 00 )
1491!
1492!--       Define x-axis (for scalar position)
1493          CALL netcdf_create_dim( id_set_3d(av), 'x', nx+1, id_dim_x_3d(av),   &
1494                                  73 )
1495          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av) /), 'x',   &
1496                                  NF90_DOUBLE, id_var_x_3d(av), 'meters', '',  &
1497                                  74, 75, 00 )
1498!
1499!--       Define x-axis (for u position)
1500          CALL netcdf_create_dim( id_set_3d(av), 'xu', nx+1, id_dim_xu_3d(av), &
1501                                  358 )
1502          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_xu_3d(av) /), 'xu', &
1503                                  NF90_DOUBLE, id_var_xu_3d(av), 'meters', '', &
1504                                  359, 360, 000 )
1505!
1506!--       Define y-axis (for scalar position)
1507          CALL netcdf_create_dim( id_set_3d(av), 'y', ny+1, id_dim_y_3d(av),   &
1508                                  76 )
1509          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_y_3d(av) /), 'y',   &
1510                                  NF90_DOUBLE, id_var_y_3d(av), 'meters', '',  &
1511                                  77, 78, 00 )
1512!
1513!--       Define y-axis (for v position)
1514          CALL netcdf_create_dim( id_set_3d(av), 'yv', ny+1, id_dim_yv_3d(av), &
1515                                  361 )
1516          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_yv_3d(av) /), 'yv', &
1517                                  NF90_DOUBLE, id_var_yv_3d(av), 'meters', '', &
1518                                  362, 363, 000 )
1519!
1520!--       Define UTM and geographic coordinates
1521          CALL define_geo_coordinates( id_set_3d(av),         &
1522                  (/ id_dim_x_3d(av), id_dim_xu_3d(av) /),    &
1523                  (/ id_dim_y_3d(av), id_dim_yv_3d(av) /),    &
1524                  id_var_eutm_3d(:,av), id_var_nutm_3d(:,av), &
1525                  id_var_lat_3d(:,av), id_var_lon_3d(:,av)    )
1526!
1527!--       Define coordinate-reference system
1528          CALL netcdf_create_crs( id_set_3d(av), 000 )
1529!
1530!--       In case of non-flat topography define 2d-arrays containing the height
1531!--       information. Only output 2d topography information in case of parallel
1532!--       output.
1533          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
1534               netcdf_data_format > 4 )  THEN
1535!
1536!--          Define zusi = zu(nzb_s_inner)
1537             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av),        &
1538                                     id_dim_y_3d(av) /), 'zusi', NF90_DOUBLE,  &
1539                                     id_var_zusi_3d(av), 'meters',             &
1540                                     'zu(nzb_s_inner)', 413, 414, 415 )
1541!
1542!--          Define zwwi = zw(nzb_w_inner)
1543             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av),        &
1544                                     id_dim_y_3d(av) /), 'zwwi', NF90_DOUBLE,  &
1545                                     id_var_zwwi_3d(av), 'meters',             &
1546                                     'zw(nzb_w_inner)', 416, 417, 418 )
1547
1548          ENDIF
1549
1550          IF ( land_surface )  THEN
1551!
1552!--          Define vertical coordinate grid (zs grid)
1553             CALL netcdf_create_dim( id_set_3d(av), 'zs_3d',                   &
1554                                     nzt_soil-nzb_soil+1, id_dim_zs_3d(av), 70 )
1555             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zs_3d(av) /),    &
1556                                     'zs_3d', NF90_DOUBLE, id_var_zs_3d(av),   &
1557                                     'meters', '', 71, 72, 00 )
1558
1559          ENDIF
1560
1561          IF ( plant_canopy )  THEN
1562!
1563!--          Define vertical coordinate grid (zpc grid)
1564             CALL netcdf_create_dim( id_set_3d(av), 'zpc_3d',                  &
1565                                     pch_index+1, id_dim_zpc_3d(av), 70 )
1566             !netcdf_create_dim(ncid, dim_name, ncdim_type, ncdim_id, error_no)
1567             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zpc_3d(av) /),    &
1568                                     'zpc_3d', NF90_DOUBLE, id_var_zpc_3d(av),   &
1569                                     'meters', '', 71, 72, 00 )
1570
1571          ENDIF
1572
1573!
1574!--       Define the variables
1575          var_list = ';'
1576          i = 1
1577
1578          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
1579!
1580!--          Temporary solution to account for data output within the new urban
1581!--          surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
1582             trimvar = TRIM( do3d(av,i) )
1583             IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
1584                trimvar = 'usm_output'
1585             ENDIF
1586!
1587!--          Check for the grid
1588             found = .FALSE.
1589             SELECT CASE ( trimvar )
1590!
1591!--             Most variables are defined on the scalar grid
1592                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',   &
1593                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
1594                       's', 'theta', 'thetal', 'thetav' )
1595
1596                   grid_x = 'x'
1597                   grid_y = 'y'
1598                   grid_z = 'zu'
1599!
1600!--             u grid
1601                CASE ( 'u' )
1602
1603                   grid_x = 'xu'
1604                   grid_y = 'y'
1605                   grid_z = 'zu'
1606!
1607!--             v grid
1608                CASE ( 'v' )
1609
1610                   grid_x = 'x'
1611                   grid_y = 'yv'
1612                   grid_z = 'zu'
1613!
1614!--             w grid
1615                CASE ( 'w' )
1616
1617                   grid_x = 'x'
1618                   grid_y = 'y'
1619                   grid_z = 'zw'
1620
1621!
1622!--             Block of urban surface model outputs
1623                CASE ( 'usm_output' )
1624                   CALL usm_define_netcdf_grid( do3d(av,i), found, &
1625                                                   grid_x, grid_y, grid_z )
1626
1627                CASE DEFAULT
1628
1629                   CALL tcm_define_netcdf_grid( do3d(av,i), found, &
1630                                                   grid_x, grid_y, grid_z )
1631
1632!
1633!--                Check for land surface quantities
1634                   IF ( .NOT. found .AND. land_surface )  THEN
1635                      CALL lsm_define_netcdf_grid( do3d(av,i), found, grid_x,  &
1636                                                   grid_y, grid_z )
1637                   ENDIF
1638!
1639!--                Check for ocean quantities
1640                   IF ( .NOT. found  .AND.  ocean_mode )  THEN
1641                      CALL ocean_define_netcdf_grid( do3d(av,i), found,  &
1642                                                     grid_x, grid_y, grid_z )
1643                   ENDIF
1644
1645!
1646!--                Check for plant canopy quantities
1647                   IF ( .NOT. found  .AND.  plant_canopy )  THEN
1648                      CALL pcm_define_netcdf_grid( do3d(av,i), found, grid_x,  &
1649                                                   grid_y, grid_z )
1650                   ENDIF
1651
1652!
1653!--                Check for radiation quantities
1654                   IF ( .NOT. found  .AND.  radiation )  THEN
1655                      CALL radiation_define_netcdf_grid( do3d(av,i), found,    &
1656                                                         grid_x, grid_y,       &
1657                                                         grid_z )
1658                   ENDIF
1659
1660!--                Check for gust module quantities
1661                   IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
1662                      CALL gust_define_netcdf_grid( do3d(av,i), found, grid_x, &
1663                                                    grid_y, grid_z )
1664                   ENDIF
1665!
1666!--                Check for indoor model quantities
1667                   IF ( .NOT. found .AND. indoor_model ) THEN
1668                      CALL im_define_netcdf_grid( do3d(av,i), found,           &
1669                                                  grid_x, grid_y, grid_z )
1670                   ENDIF
1671
1672!
1673!--                Check for biometeorology quantities
1674                   IF ( .NOT. found  .AND.  biometeorology )  THEN
1675                      CALL bio_define_netcdf_grid( do3d(av,i), found,          &
1676                                                   grid_x, grid_y, grid_z )
1677                   ENDIF
1678
1679!
1680!--                Check for chemistry quantities
1681                   IF ( .NOT. found  .AND.  air_chemistry )  THEN
1682                      CALL chem_define_netcdf_grid( do3d(av,i), found,         &
1683                                                    grid_x, grid_y, grid_z )
1684                   ENDIF
1685
1686!
1687!--                Check for SALSA quantities
1688                   IF ( .NOT. found  .AND.  salsa )  THEN
1689                      CALL salsa_define_netcdf_grid( do3d(av,i), found, grid_x,&
1690                                                     grid_y, grid_z )
1691                   ENDIF
1692!
1693!--                Check for user-defined quantities
1694                   IF ( .NOT. found  .AND.  user_module_enabled )  THEN
1695                      CALL user_define_netcdf_grid( do3d(av,i), found, grid_x, &
1696                                                    grid_y, grid_z )
1697                   ENDIF
1698
1699                   IF ( .NOT. found )                                          &
1700                      CALL doq_define_netcdf_grid( do3d(av,i), found, grid_x,  &
1701                                                   grid_y, grid_z        )
1702
1703                   IF ( .NOT. found )  THEN
1704                      WRITE ( message_string, * ) 'no grid defined for varia', &
1705                                                  'ble ', TRIM( do3d(av,i) )
1706                      CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, &
1707                                    6, 0 )
1708                   ENDIF
1709
1710             END SELECT
1711
1712!
1713!--          Select the respective dimension ids
1714             IF ( grid_x == 'x' )  THEN
1715                id_x = id_dim_x_3d(av)
1716             ELSEIF ( grid_x == 'xu' )  THEN
1717                id_x = id_dim_xu_3d(av)
1718             ENDIF
1719
1720             IF ( grid_y == 'y' )  THEN
1721                id_y = id_dim_y_3d(av)
1722             ELSEIF ( grid_y == 'yv' )  THEN
1723                id_y = id_dim_yv_3d(av)
1724             ENDIF
1725
1726             IF ( grid_z == 'zu' )  THEN
1727                id_z = id_dim_zu_3d(av)
1728             ELSEIF ( grid_z == 'zw' )  THEN
1729                id_z = id_dim_zw_3d(av)
1730             ELSEIF ( grid_z == 'zs' )  THEN
1731                id_z = id_dim_zs_3d(av)
1732             ELSEIF ( grid_z == 'zpc' )  THEN
1733                id_z = id_dim_zpc_3d(av)
1734             ENDIF
1735
1736!
1737!--          Define the grid
1738             CALL netcdf_create_var( id_set_3d(av),(/ id_x, id_y, id_z,        &
1739                                     id_dim_time_3d(av) /), do3d(av,i),        &
1740                                     nc_precision(4), id_var_do3d(av,i),       &
1741                                     TRIM( do3d_unit(av,i) ), do3d(av,i), 79,  &
1742                                     80, 357, .TRUE. )
1743#if defined( __netcdf4_parallel )
1744             IF ( netcdf_data_format > 4 )  THEN
1745!
1746!--             Set no fill for every variable to increase performance.
1747                nc_stat = NF90_DEF_VAR_FILL( id_set_3d(av),     &
1748                                             id_var_do3d(av,i), &
1749                                             NF90_NOFILL, 0 )
1750                CALL netcdf_handle_error( 'netcdf_define_header', 532 )
1751!
1752!--             Set collective io operations for parallel io
1753                nc_stat = NF90_VAR_PAR_ACCESS( id_set_3d(av),     &
1754                                               id_var_do3d(av,i), &
1755                                               NF90_COLLECTIVE )
1756                CALL netcdf_handle_error( 'netcdf_define_header', 445 )
1757             ENDIF
1758#endif
1759             var_list = TRIM( var_list ) // TRIM( do3d(av,i) ) // ';'
1760
1761             i = i + 1
1762
1763          ENDDO
1764
1765!
1766!--       No arrays to output
1767          IF ( i == 1 )  RETURN
1768
1769!
1770!--       Write the list of variables as global attribute (this is used by
1771!--       restart runs and by combine_plot_fields)
1772          nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
1773                                  var_list )
1774          CALL netcdf_handle_error( 'netcdf_define_header', 81 )
1775
1776!
1777!--       Set general no fill, otherwise the performance drops significantly for
1778!--       parallel output.
1779          nc_stat = NF90_SET_FILL( id_set_3d(av), NF90_NOFILL, oldmode )
1780          CALL netcdf_handle_error( 'netcdf_define_header', 528 )
1781
1782!
1783!--       Leave netCDF define mode
1784          nc_stat = NF90_ENDDEF( id_set_3d(av) )
1785          CALL netcdf_handle_error( 'netcdf_define_header', 82 )
1786
1787!
1788!--       These data are only written by PE0 for parallel output to increase
1789!--       the performance.
1790          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
1791!
1792!--          Write data for x (shifted by +dx/2) and xu axis
1793             ALLOCATE( netcdf_data(0:nx) )
1794
1795             DO  i = 0, nx
1796                netcdf_data(i) = ( i + 0.5 ) * dx
1797             ENDDO
1798
1799             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_x_3d(av),  &
1800                                     netcdf_data, start = (/ 1 /),    &
1801                                     count = (/ nx+1 /) )
1802             CALL netcdf_handle_error( 'netcdf_define_header', 83 )
1803
1804             DO  i = 0, nx
1805                netcdf_data(i) = i * dx
1806             ENDDO
1807
1808             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_xu_3d(av), &
1809                                     netcdf_data, start = (/ 1 /),    &
1810                                     count = (/ nx+1 /) )
1811             CALL netcdf_handle_error( 'netcdf_define_header', 385 )
1812
1813             DEALLOCATE( netcdf_data )
1814
1815!
1816!--          Write data for y (shifted by +dy/2) and yv axis
1817             ALLOCATE( netcdf_data(0:ny) )
1818
1819             DO  i = 0, ny
1820                netcdf_data(i) = ( i + 0.5_wp ) * dy
1821             ENDDO
1822
1823             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_y_3d(av),  &
1824                                     netcdf_data, start = (/ 1 /),    &
1825                                     count = (/ ny+1 /) )
1826             CALL netcdf_handle_error( 'netcdf_define_header', 84 )
1827
1828             DO  i = 0, ny
1829                netcdf_data(i) = i * dy
1830             ENDDO
1831
1832             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_yv_3d(av), &
1833                                     netcdf_data, start = (/ 1 /),    &
1834                                     count = (/ ny+1 /))
1835             CALL netcdf_handle_error( 'netcdf_define_header', 387 )
1836
1837             DEALLOCATE( netcdf_data )
1838
1839!
1840!--          Write UTM coordinates
1841             IF ( init_model%rotation_angle == 0.0_wp )  THEN
1842!
1843!--             1D in case of no rotation
1844                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
1845!
1846!--             x coordinates
1847                ALLOCATE( netcdf_data(0:nx) )
1848                DO  k = 0, 2
1849!
1850!--                Scalar grid points
1851                   IF ( k == 0 )  THEN
1852                      shift_x = 0.5
1853!
1854!--                u grid points
1855                   ELSEIF ( k == 1 )  THEN
1856                      shift_x = 0.0
1857!
1858!--                v grid points
1859                   ELSEIF ( k == 2 )  THEN
1860                      shift_x = 0.5
1861                   ENDIF
1862
1863                   DO  i = 0, nx
1864                     netcdf_data(i) = init_model%origin_x            &
1865                                    + cos_ra * ( i + shift_x ) * dx
1866                   ENDDO
1867
1868                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_eutm_3d(k,av),&
1869                                           netcdf_data, start = (/ 1 /),   &
1870                                           count = (/ nx+1 /) )
1871                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
1872
1873                ENDDO
1874                DEALLOCATE( netcdf_data )
1875!
1876!--             y coordinates
1877                ALLOCATE( netcdf_data(0:ny) )
1878                DO  k = 0, 2
1879!
1880!--                Scalar grid points
1881                   IF ( k == 0 )  THEN
1882                      shift_y = 0.5
1883!
1884!--                u grid points
1885                   ELSEIF ( k == 1 )  THEN
1886                      shift_y = 0.5
1887!
1888!--                v grid points
1889                   ELSEIF ( k == 2 )  THEN
1890                      shift_y = 0.0
1891                   ENDIF
1892
1893                   DO  j = 0, ny
1894                      netcdf_data(j) = init_model%origin_y            &
1895                                     + cos_ra * ( j + shift_y ) * dy
1896                   ENDDO
1897
1898                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_nutm_3d(k,av),&
1899                                           netcdf_data, start = (/ 1 /),   &
1900                                           count = (/ ny+1 /) )
1901                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
1902
1903                ENDDO
1904                DEALLOCATE( netcdf_data )
1905
1906             ELSE
1907!
1908!--             2D in case of rotation
1909                ALLOCATE( netcdf_data_2d(0:nx,0:ny) )
1910                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
1911                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
1912
1913                DO  k = 0, 2
1914!
1915!--               Scalar grid points
1916                  IF ( k == 0 )  THEN
1917                     shift_x = 0.5 ; shift_y = 0.5
1918!
1919!--               u grid points
1920                  ELSEIF ( k == 1 )  THEN
1921                     shift_x = 0.0 ; shift_y = 0.5
1922!
1923!--               v grid points
1924                  ELSEIF ( k == 2 )  THEN
1925                     shift_x = 0.5 ; shift_y = 0.0
1926                  ENDIF
1927
1928                  DO  j = 0, ny
1929                     DO  i = 0, nx
1930                        netcdf_data_2d(i,j) = init_model%origin_x            &
1931                                            + cos_ra * ( i + shift_x ) * dx  &
1932                                            + sin_ra * ( j + shift_y ) * dy
1933                     ENDDO
1934                  ENDDO
1935
1936                  nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_eutm_3d(k,av),  &
1937                                          netcdf_data_2d, start = (/ 1, 1 /),   &
1938                                          count = (/ nx+1, ny+1 /) )
1939                  CALL netcdf_handle_error( 'netcdf_define_header', 555 )
1940
1941                  DO  j = 0, ny
1942                     DO  i = 0, nx
1943                        netcdf_data_2d(i,j) = init_model%origin_y            &
1944                                            - sin_ra * ( i + shift_x ) * dx  &
1945                                            + cos_ra * ( j + shift_y ) * dy
1946                     ENDDO
1947                  ENDDO
1948
1949                  nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_nutm_3d(k,av),  &
1950                                          netcdf_data_2d, start = (/ 1, 1 /),   &
1951                                          count = (/ nx+1, ny+1 /) )
1952                  CALL netcdf_handle_error( 'netcdf_define_header', 556 )
1953
1954                ENDDO
1955                DEALLOCATE( netcdf_data_2d )
1956             ENDIF
1957!
1958!--          Write zu and zw data (vertical axes)
1959             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zu_3d(av),  &
1960                                     zu(nzb:nz_do3d), start = (/ 1 /), &
1961                                     count = (/ nz_do3d-nzb+1 /) )
1962             CALL netcdf_handle_error( 'netcdf_define_header', 85 )
1963
1964
1965             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zw_3d(av),  &
1966                                     zw(nzb:nz_do3d), start = (/ 1 /), &
1967                                     count = (/ nz_do3d-nzb+1 /) )
1968             CALL netcdf_handle_error( 'netcdf_define_header', 86 )
1969
1970             IF ( land_surface )  THEN
1971!
1972!--             Write zs grid
1973                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zs_3d(av),  &
1974                                        - zs(nzb_soil:nzt_soil), start = (/ 1 /), &
1975                                        count = (/ nzt_soil-nzb_soil+1 /) )
1976                CALL netcdf_handle_error( 'netcdf_define_header', 86 )
1977             ENDIF
1978
1979             IF ( plant_canopy )  THEN
1980!
1981!--             Write zpc grid
1982                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zpc_3d(av),  &
1983                                        zu(nzb:nzb+pch_index), start = (/ 1 /), &
1984                                        count = (/ pch_index+1 /) )
1985                CALL netcdf_handle_error( 'netcdf_define_header', 86 )
1986             ENDIF
1987
1988          ENDIF
1989!
1990!--       Write lon and lat data. Only for parallel output.
1991          IF ( netcdf_data_format > 4 )  THEN
1992
1993             ALLOCATE( lat(nxl:nxr,nys:nyn) )
1994             ALLOCATE( lon(nxl:nxr,nys:nyn) )
1995             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
1996             sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
1997
1998             DO  k = 0, 2
1999!
2000!--             Scalar grid points
2001                IF ( k == 0 )  THEN
2002                   shift_x = 0.5 ; shift_y = 0.5
2003!
2004!--             u grid points
2005                ELSEIF ( k == 1 )  THEN
2006                   shift_x = 0.0 ; shift_y = 0.5
2007!
2008!--             v grid points
2009                ELSEIF ( k == 2 )  THEN
2010                   shift_x = 0.5 ; shift_y = 0.0
2011                ENDIF
2012
2013                DO  j = nys, nyn
2014                   DO  i = nxl, nxr
2015                      eutm = init_model%origin_x            &
2016                           + cos_ra * ( i + shift_x ) * dx  &
2017                           + sin_ra * ( j + shift_y ) * dy
2018                      nutm = init_model%origin_y            &
2019                           - sin_ra * ( i + shift_x ) * dx  &
2020                           + cos_ra * ( j + shift_y ) * dy
2021
2022                      CALL  convert_utm_to_geographic( crs_list,          &
2023                                                       eutm, nutm,        &
2024                                                       lon(i,j), lat(i,j) )
2025                   ENDDO
2026                ENDDO
2027
2028                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_lon_3d(k,av), &
2029                                     lon, start = (/ nxl+1, nys+1 /),       &
2030                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
2031                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
2032
2033                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_lat_3d(k,av), &
2034                                     lat, start = (/ nxl+1, nys+1 /),       &
2035                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
2036                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
2037             ENDDO
2038
2039             DEALLOCATE( lat )
2040             DEALLOCATE( lon )
2041
2042          ENDIF
2043!
2044!--       In case of non-flat topography write height information. Only for
2045!--       parallel netcdf output.
2046          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
2047               netcdf_data_format > 4 )  THEN
2048
2049!             IF ( nxr == nx  .AND.  nyn /= ny )  THEN
2050!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av),     &
2051!                                        zu_s_inner(nxl:nxr+1,nys:nyn),         &
2052!                                        start = (/ nxl+1, nys+1 /),            &
2053!                                        count = (/ nxr-nxl+2, nyn-nys+1 /) )
2054!             ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
2055!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av),     &
2056!                                        zu_s_inner(nxl:nxr,nys:nyn+1),         &
2057!                                        start = (/ nxl+1, nys+1 /),            &
2058!                                        count = (/ nxr-nxl+1, nyn-nys+2 /) )
2059!             ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
2060!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av),     &
2061!                                        zu_s_inner(nxl:nxr+1,nys:nyn+1),       &
2062!                                        start = (/ nxl+1, nys+1 /),            &
2063!                                        count = (/ nxr-nxl+2, nyn-nys+2 /) )
2064!             ELSE
2065                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av),     &
2066                                        zu_s_inner(nxl:nxr,nys:nyn),           &
2067                                        start = (/ nxl+1, nys+1 /),            &
2068                                        count = (/ nxr-nxl+1, nyn-nys+1 /) )
2069!             ENDIF
2070             CALL netcdf_handle_error( 'netcdf_define_header', 419 )
2071
2072!             IF ( nxr == nx  .AND.  nyn /= ny )  THEN
2073!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av),     &
2074!                                        zw_w_inner(nxl:nxr+1,nys:nyn),         &
2075!                                        start = (/ nxl+1, nys+1 /),            &
2076!                                        count = (/ nxr-nxl+2, nyn-nys+1 /) )
2077!             ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
2078!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av),     &
2079!                                        zw_w_inner(nxl:nxr,nys:nyn+1),         &
2080!                                        start = (/ nxl+1, nys+1 /),            &
2081!                                        count = (/ nxr-nxl+1, nyn-nys+2 /) )
2082!             ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
2083!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av),     &
2084!                                        zw_w_inner(nxl:nxr+1,nys:nyn+1),       &
2085!                                        start = (/ nxl+1, nys+1 /),            &
2086!                                        count = (/ nxr-nxl+2, nyn-nys+2 /) )
2087!             ELSE
2088                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av),     &
2089                                        zw_w_inner(nxl:nxr,nys:nyn),           &
2090                                        start = (/ nxl+1, nys+1 /),            &
2091                                        count = (/ nxr-nxl+1, nyn-nys+1 /) )
2092!             ENDIF
2093             CALL netcdf_handle_error( 'netcdf_define_header', 420 )
2094
2095          ENDIF
2096
2097       CASE ( '3d_ext' )
2098
2099!
2100!--       Get the list of variables and compare with the actual run.
2101!--       First var_list_old has to be reset, since GET_ATT does not assign
2102!--       trailing blanks.
2103          var_list_old = ' '
2104          nc_stat = NF90_GET_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
2105                                  var_list_old )
2106          CALL netcdf_handle_error( 'netcdf_define_header', 87 )
2107
2108          var_list = ';'
2109          i = 1
2110          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
2111             var_list = TRIM(var_list) // TRIM( do3d(av,i) ) // ';'
2112             i = i + 1
2113          ENDDO
2114
2115          IF ( av == 0 )  THEN
2116             var = '(3d)'
2117          ELSE
2118             var = '(3d_av)'
2119          ENDIF
2120
2121          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
2122             message_string = 'netCDF file for volume data ' //             &
2123                              TRIM( var ) // ' from previous run found,' // &
2124                              '&but this file cannot be extended due to' // &
2125                              ' variable mismatch.' //                      &
2126                              '&New file is created instead.'
2127             CALL message( 'define_netcdf_header', 'PA0245', 0, 1, 0, 6, 0 )
2128             extend = .FALSE.
2129             RETURN
2130          ENDIF
2131
2132!
2133!--       Get and compare the number of vertical gridpoints
2134          nc_stat = NF90_INQ_VARID( id_set_3d(av), 'zu_3d', id_var_zu_3d(av) )
2135          CALL netcdf_handle_error( 'netcdf_define_header', 88 )
2136
2137          nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_zu_3d(av), &
2138                                           dimids = id_dim_zu_3d_old )
2139          CALL netcdf_handle_error( 'netcdf_define_header', 89 )
2140          id_dim_zu_3d(av) = id_dim_zu_3d_old(1)
2141
2142          nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_zu_3d(av), &
2143                                            len = nz_old )
2144          CALL netcdf_handle_error( 'netcdf_define_header', 90 )
2145
2146          IF ( nz_do3d-nzb+1 /= nz_old )  THEN
2147              message_string = 'netCDF file for volume data ' //             &
2148                               TRIM( var ) // ' from previous run found,' // &
2149                               '&but this file cannot be extended due to' // &
2150                               ' mismatch in number of' //                   &
2151                               ' vertical grid points (nz_do3d).' //         &
2152                               '&New file is created instead.'
2153             CALL message( 'define_netcdf_header', 'PA0246', 0, 1, 0, 6, 0 )
2154             extend = .FALSE.
2155             RETURN
2156          ENDIF
2157
2158!
2159!--       Get the id of the time coordinate (unlimited coordinate) and its
2160!--       last index on the file. The next time level is pl3d..count+1.
2161!--       The current time must be larger than the last output time
2162!--       on the file.
2163          nc_stat = NF90_INQ_VARID( id_set_3d(av), 'time', id_var_time_3d(av) )
2164          CALL netcdf_handle_error( 'netcdf_define_header', 91 )
2165
2166          nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_time_3d(av), &
2167                                           dimids = id_dim_time_old )
2168          CALL netcdf_handle_error( 'netcdf_define_header', 92 )
2169
2170          id_dim_time_3d(av) = id_dim_time_old(1)
2171
2172          nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_time_3d(av), &
2173                                            len = ntime_count )
2174          CALL netcdf_handle_error( 'netcdf_define_header', 93 )
2175
2176!
2177!--       For non-parallel output use the last output time level of the netcdf
2178!--       file because the time dimension is unlimited. In case of parallel
2179!--       output the variable ntime_count could get the value of 9*10E36 because
2180!--       the time dimension is limited.
2181          IF ( netcdf_data_format < 5 ) do3d_time_count(av) = ntime_count
2182
2183          nc_stat = NF90_GET_VAR( id_set_3d(av), id_var_time_3d(av), &
2184                                  last_time_coordinate,              &
2185                                  start = (/ do3d_time_count(av) /), &
2186                                  count = (/ 1 /) )
2187          CALL netcdf_handle_error( 'netcdf_define_header', 94 )
2188
2189          IF ( last_time_coordinate(1) >= simulated_time )  THEN
2190             message_string = 'netCDF file for volume data ' //             &
2191                              TRIM( var ) // ' from previous run found,' // &
2192                              '&but this file cannot be extended becaus' // &
2193                              'e the current output time' //                &
2194                              '&is less or equal than the last output t' // &
2195                              'ime on this file.' //                        &
2196                              '&New file is created instead.'
2197             CALL message( 'define_netcdf_header', 'PA0247', 0, 1, 0, 6, 0 )
2198             do3d_time_count(av) = 0
2199             extend = .FALSE.
2200             RETURN
2201          ENDIF
2202
2203          IF ( netcdf_data_format > 4 )  THEN
2204!
2205!--          Check if the needed number of output time levels is increased
2206!--          compared to the number of time levels in the existing file.
2207             IF ( ntdim_3d(av) > ntime_count )  THEN
2208                message_string = 'netCDF file for volume data ' // &
2209                                 TRIM( var ) // ' from previous run found,' // &
2210                                 '&but this file cannot be extended becaus' // &
2211                                 'e the number of output time levels has b' // &
2212                                 'een increased compared to the previous s' // &
2213                                 'imulation.' //                               &
2214                                 '&New file is created instead.'
2215                CALL message( 'define_netcdf_header', 'PA0388', 0, 1, 0, 6, 0 )
2216                do3d_time_count(av) = 0
2217                extend = .FALSE.
2218!
2219!--             Recalculate the needed time levels for the new file.
2220                IF ( av == 0 )  THEN
2221                   ntdim_3d(0) = CEILING(                               &
2222                           ( end_time - MAX( skip_time_do3d,            &
2223                                             simulated_time_at_begin )  &
2224                           ) / dt_do3d )
2225                   IF ( do3d_at_begin )  ntdim_3d(0) = ntdim_3d(0) + 1
2226                ELSE
2227                   ntdim_3d(1) = CEILING(                               &
2228                           ( end_time - MAX( skip_time_data_output_av,  &
2229                                             simulated_time_at_begin )  &
2230                           ) / dt_data_output_av )
2231                ENDIF
2232                RETURN
2233             ENDIF
2234          ENDIF
2235
2236!
2237!--       Dataset seems to be extendable.
2238!--       Now get the variable ids.
2239          i = 1
2240          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
2241             nc_stat = NF90_INQ_VARID( id_set_3d(av), TRIM( do3d(av,i) ), &
2242                                       id_var_do3d(av,i) )
2243             CALL netcdf_handle_error( 'netcdf_define_header', 95 )
2244#if defined( __netcdf4_parallel )
2245!
2246!--          Set collective io operations for parallel io
2247             IF ( netcdf_data_format > 4 )  THEN
2248                nc_stat = NF90_VAR_PAR_ACCESS( id_set_3d(av),     &
2249                                               id_var_do3d(av,i), &
2250                                               NF90_COLLECTIVE )
2251                CALL netcdf_handle_error( 'netcdf_define_header', 453 )
2252             ENDIF
2253#endif
2254             i = i + 1
2255          ENDDO
2256
2257!
2258!--       Update the title attribute on file
2259!--       In order to avoid 'data mode' errors if updated attributes are larger
2260!--       than their original size, NF90_PUT_ATT is called in 'define mode'
2261!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
2262!--       performance loss due to data copying; an alternative strategy would be
2263!--       to ensure equal attribute size. Maybe revise later.
2264          IF ( av == 0 )  THEN
2265             time_average_text = ' '
2266          ELSE
2267             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
2268                                                            averaging_interval
2269          ENDIF
2270          nc_stat = NF90_REDEF( id_set_3d(av) )
2271          CALL netcdf_handle_error( 'netcdf_define_header', 429 )
2272          nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'title', &
2273                                  TRIM( run_description_header ) //    &
2274                                  TRIM( time_average_text ) )
2275          CALL netcdf_handle_error( 'netcdf_define_header', 96 )
2276          nc_stat = NF90_ENDDEF( id_set_3d(av) )
2277          CALL netcdf_handle_error( 'netcdf_define_header', 430 )
2278          message_string = 'netCDF file for volume data ' //             &
2279                           TRIM( var ) // ' from previous run found.' // &
2280                           '&This file will be extended.'
2281          CALL message( 'define_netcdf_header', 'PA0248', 0, 0, 0, 6, 0 )
2282
2283
2284       CASE ( 'ag_new' )
2285
2286!
2287!--       Define some global attributes of the dataset
2288          nc_stat = NF90_PUT_ATT( id_set_agt, NF90_GLOBAL, 'title', &
2289                                  TRIM( run_description_header ) )
2290          CALL netcdf_handle_error( 'netcdf_define_header', 330 )
2291!
2292!--       Switch for unlimited time dimension
2293          IF ( agent_time_unlimited ) THEN
2294             CALL netcdf_create_dim( id_set_agt, 'time', NF90_UNLIMITED,       &
2295                                     id_dim_time_agt, 331 )
2296          ELSE
2297             CALL netcdf_create_dim( id_set_agt, 'time',                       &
2298                                     INT( ( MIN( multi_agent_system_end,       &
2299                                                 end_time ) -                  &
2300                                            multi_agent_system_start ) /       &
2301                                            dt_write_agent_data * 1.1 ),       &
2302                                     id_dim_time_agt, 331 )
2303          ENDIF
2304
2305          CALL netcdf_create_var( id_set_agt, (/ id_dim_time_agt /), 'time',   &
2306                                  NF90_REAL4, id_var_time_agt, 'seconds', 'time',  &
2307                                  332, 333, 000 )
2308          CALL netcdf_create_att( id_set_agt, id_var_time_agt, 'standard_name', 'time', 000)
2309          CALL netcdf_create_att( id_set_agt, id_var_time_agt, 'axis', 'T', 000)
2310
2311          CALL netcdf_create_dim( id_set_agt, 'agent_number',                  &
2312                                  dim_size_agtnum, id_dim_agtnum, 334 )
2313
2314          CALL netcdf_create_var( id_set_agt, (/ id_dim_agtnum /),             &
2315                                  'agent_number', NF90_REAL4,                  &
2316                                  id_var_agtnum, 'agent number', '', 335,      &
2317                                  336, 000 )
2318!
2319!--       Define variable which contains the real number of agents in use
2320          CALL netcdf_create_var( id_set_agt, (/ id_dim_time_agt /),           &
2321                                  'real_num_of_agt', NF90_REAL4,               &
2322                                  id_var_rnoa_agt, 'agent number', '', 337,    &
2323                                  338, 000 )
2324          i = 1
2325          CALL netcdf_create_var( id_set_agt, (/ id_dim_agtnum,                &
2326                                  id_dim_time_agt /), agt_var_names(i),        &
2327                                  NF90_DOUBLE, id_var_agt(i),                  &
2328                                  TRIM( agt_var_units(i) ),                    &
2329                                  TRIM( agt_var_names(i) ), 339, 340, 341 )
2330!
2331!--       Define the variables
2332          DO  i = 2, 6
2333             CALL netcdf_create_var( id_set_agt, (/ id_dim_agtnum,             &
2334                                     id_dim_time_agt /), agt_var_names(i),     &
2335                                     NF90_REAL4, id_var_agt(i),                &
2336                                     TRIM( agt_var_units(i) ),                 &
2337                                     TRIM( agt_var_names(i) ), 339, 340, 341 )
2338
2339          ENDDO
2340!
2341!--       Define vars for biometeorology
2342          CALL netcdf_create_var( id_set_agt, (/ id_dim_agtnum,                &
2343                                  id_dim_time_agt /), agt_var_names(9),        &
2344                                  nc_precision(8), id_var_agt(9),              &
2345                                  TRIM( agt_var_units(9) ),                    &
2346                                  TRIM( agt_var_names(9) ), 339, 340, 341 )
2347
2348!
2349!--       Leave netCDF define mode
2350          nc_stat = NF90_ENDDEF( id_set_agt )
2351          CALL netcdf_handle_error( 'netcdf_define_header', 342 )
2352
2353
2354!        CASE ( 'ag_ext' )
2355! !+?agent extend output for restart runs has to be adapted
2356!
2357! !
2358! !--       Get the id of the time coordinate (unlimited coordinate) and its
2359! !--       last index on the file. The next time level is prt..count+1.
2360! !--       The current time must be larger than the last output time
2361! !--       on the file.
2362!           nc_stat = NF90_INQ_VARID( id_set_agt, 'time', id_var_time_agt )
2363!           CALL netcdf_handle_error( 'netcdf_define_header', 343 )
2364!
2365!           nc_stat = NF90_INQUIRE_VARIABLE( id_set_agt, id_var_time_agt, &
2366!                                            dimids = id_dim_time_old )
2367!           CALL netcdf_handle_error( 'netcdf_define_header', 344 )
2368!           id_dim_time_agt = id_dim_time_old(1)
2369!
2370!           nc_stat = NF90_INQUIRE_DIMENSION( id_set_agt, id_dim_time_agt, &
2371!                                             len = agt_time_count )
2372!           CALL netcdf_handle_error( 'netcdf_define_header', 345 )
2373!
2374!           nc_stat = NF90_GET_VAR( id_set_agt, id_var_time_agt,  &
2375!                                   last_time_coordinate,         &
2376!                                   start = (/ agt_time_count /), &
2377!                                   count = (/ 1 /) )
2378!           CALL netcdf_handle_error( 'netcdf_define_header', 346 )
2379!
2380!           IF ( last_time_coordinate(1) >= simulated_time )  THEN
2381!              message_string = 'netCDF file for agents ' //                  &
2382!                               'from previous run found,' //                 &
2383!                               '&but this file cannot be extended becaus' // &
2384!                               'e the current output time' //                &
2385!                               '&is less or equal than the last output t' // &
2386!                               'ime on this file.' //                        &
2387!                               '&New file is created instead.'
2388!              CALL message( 'define_netcdf_header', 'PA0265', 0, 1, 0, 6, 0 )
2389!              agt_time_count = 0
2390!              extend = .FALSE.
2391!              RETURN
2392!           ENDIF
2393!
2394! !
2395! !--       Dataset seems to be extendable.
2396! !--       Now get the variable ids.
2397!           nc_stat = NF90_INQ_VARID( id_set_agt, 'real_num_of_agt', &
2398!                                     id_var_rnoa_agt )
2399!           CALL netcdf_handle_error( 'netcdf_define_header', 347 )
2400!
2401!           DO  i = 1, 17
2402!
2403!              nc_stat = NF90_INQ_VARID( id_set_agt, agt_var_names(i), &
2404!                                        id_var_prt(i) )
2405!              CALL netcdf_handle_error( 'netcdf_define_header', 348 )
2406!
2407!           ENDDO
2408!
2409!           message_string = 'netCDF file for particles ' // &
2410!                            'from previous run found.' //   &
2411!                            '&This file will be extended.'
2412!           CALL message( 'define_netcdf_header', 'PA0266', 0, 0, 0, 6, 0 )
2413
2414
2415       CASE ( 'xy_new' )
2416
2417!
2418!--       Define some global attributes of the dataset
2419          IF ( av == 0 )  THEN
2420             CALL netcdf_create_global_atts( id_set_xy(av), 'xy', TRIM( run_description_header ), 97 )
2421             time_average_text = ' '
2422          ELSE
2423             CALL netcdf_create_global_atts( id_set_xy(av), 'xy_av', TRIM( run_description_header ), 97 )
2424             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
2425             nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'time_avg',   &
2426                                     TRIM( time_average_text ) )
2427             CALL netcdf_handle_error( 'netcdf_define_header', 98 )
2428          ENDIF
2429
2430!
2431!--       Define time coordinate for xy sections.
2432!--       For parallel output the time dimensions has to be limited, otherwise
2433!--       the performance drops significantly.
2434          IF ( netcdf_data_format < 5 )  THEN
2435             CALL netcdf_create_dim( id_set_xy(av), 'time', NF90_UNLIMITED,    &
2436                                     id_dim_time_xy(av), 99 )
2437          ELSE
2438             CALL netcdf_create_dim( id_set_xy(av), 'time', ntdim_2d_xy(av),   &
2439                                     id_dim_time_xy(av), 524 )
2440          ENDIF
2441
2442          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_time_xy(av) /),     &
2443                                  'time', NF90_DOUBLE, id_var_time_xy(av),     &
2444                                  'seconds', 'time', 100, 101, 000 )
2445          CALL netcdf_create_att( id_set_xy(av), id_var_time_xy(av), 'standard_name', 'time', 000)
2446          CALL netcdf_create_att( id_set_xy(av), id_var_time_xy(av), 'axis', 'T', 000)
2447!
2448!--       Define the spatial dimensions and coordinates for xy-sections.
2449!--       First, determine the number of horizontal sections to be written.
2450          IF ( section(1,1) == -9999 )  THEN
2451             RETURN
2452          ELSE
2453             ns = 1
2454             DO WHILE ( section(ns,1) /= -9999  .AND.  ns <= 100 )
2455                ns = ns + 1
2456             ENDDO
2457             ns = ns - 1
2458          ENDIF
2459
2460!
2461!--       Define vertical coordinate grid (zu grid)
2462          CALL netcdf_create_dim( id_set_xy(av), 'zu_xy', ns,                  &
2463                                  id_dim_zu_xy(av), 102 )
2464          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu_xy(av) /),       &
2465                                  'zu_xy', NF90_DOUBLE, id_var_zu_xy(av),      &
2466                                  'meters', '', 103, 104, 000 )
2467!
2468!--       Define vertical coordinate grid (zw grid)
2469          CALL netcdf_create_dim( id_set_xy(av), 'zw_xy', ns,                  &
2470                                  id_dim_zw_xy(av), 105 )
2471          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zw_xy(av) /),       &
2472                                  'zw_xy', NF90_DOUBLE, id_var_zw_xy(av),      &
2473                                  'meters', '', 106, 107, 000 )
2474
2475          IF ( land_surface )  THEN
2476
2477             ns_do = 1
2478             DO WHILE ( section(ns_do,1) /= -9999  .AND.  ns_do < nzs )
2479                ns_do = ns_do + 1
2480             ENDDO
2481!
2482!--          Define vertical coordinate grid (zs grid)
2483             CALL netcdf_create_dim( id_set_xy(av), 'zs_xy', ns_do,            &
2484                                     id_dim_zs_xy(av), 539 )
2485             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zs_xy(av) /),    &
2486                                     'zs_xy', NF90_DOUBLE, id_var_zs_xy(av),   &
2487                                     'meters', '', 540, 541, 000 )
2488
2489          ENDIF
2490
2491!
2492!--       Define a pseudo vertical coordinate grid for the surface variables
2493!--       u* and t* to store their height level
2494          CALL netcdf_create_dim( id_set_xy(av), 'zu1_xy', 1,                  &
2495                                  id_dim_zu1_xy(av), 108 )
2496          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu1_xy(av) /),      &
2497                                  'zu1_xy', NF90_DOUBLE, id_var_zu1_xy(av),    &
2498                                  'meters', '', 109, 110, 000 )
2499!
2500!--       Define a variable to store the layer indices of the horizontal cross
2501!--       sections, too
2502          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu_xy(av) /),       &
2503                                  'ind_z_xy', NF90_DOUBLE,                     &
2504                                  id_var_ind_z_xy(av), 'gridpoints', '', 111,  &
2505                                  112, 000 )
2506!
2507!--       Define x-axis (for scalar position)
2508          CALL netcdf_create_dim( id_set_xy(av), 'x', nx+1, id_dim_x_xy(av),   &
2509                                  113 )
2510          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av) /), 'x',   &
2511                                  NF90_DOUBLE, id_var_x_xy(av), 'meters', '',  &
2512                                  114, 115, 000 )
2513!
2514!--       Define x-axis (for u position)
2515          CALL netcdf_create_dim( id_set_xy(av), 'xu', nx+1,                   &
2516                                  id_dim_xu_xy(av), 388 )
2517          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_xu_xy(av) /), 'xu', &
2518                                  NF90_DOUBLE, id_var_xu_xy(av), 'meters', '', &
2519                                  389, 390, 000 )
2520!
2521!--       Define y-axis (for scalar position)
2522          CALL netcdf_create_dim( id_set_xy(av), 'y', ny+1, id_dim_y_xy(av),   &
2523                                  116 )
2524          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_y_xy(av) /), 'y',   &
2525                                  NF90_DOUBLE, id_var_y_xy(av), 'meters', '',  &
2526                                  117, 118, 000 )
2527!
2528!--       Define y-axis (for scalar position)
2529          CALL netcdf_create_dim( id_set_xy(av), 'yv', ny+1,                   &
2530                                  id_dim_yv_xy(av), 364 )
2531          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_yv_xy(av) /), 'yv', &
2532                                  NF90_DOUBLE, id_var_yv_xy(av), 'meters', '', &
2533                                  365, 366, 000 )
2534!
2535!--       Define UTM and geographic coordinates
2536          CALL define_geo_coordinates( id_set_xy(av),         &
2537                  (/ id_dim_x_xy(av), id_dim_xu_xy(av) /),    &
2538                  (/ id_dim_y_xy(av), id_dim_yv_xy(av) /),    &
2539                  id_var_eutm_xy(:,av), id_var_nutm_xy(:,av), &
2540                  id_var_lat_xy(:,av), id_var_lon_xy(:,av)    )
2541!
2542!--       Define coordinate-reference system
2543          CALL netcdf_create_crs( id_set_xy(av), 000 )
2544!
2545!--       In case of non-flat topography define 2d-arrays containing the height
2546!--       information. Only for parallel netcdf output.
2547          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
2548               netcdf_data_format > 4  )  THEN
2549!
2550!--          Define zusi = zu(nzb_s_inner)
2551             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),        &
2552                                     id_dim_y_xy(av) /), 'zusi', NF90_DOUBLE,  &
2553                                     id_var_zusi_xy(av), 'meters',             &
2554                                     'zu(nzb_s_inner)', 421, 422, 423 )
2555!
2556!--          Define zwwi = zw(nzb_w_inner)
2557             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),        &
2558                                     id_dim_y_xy(av) /), 'zwwi', NF90_DOUBLE,  &
2559                                     id_var_zwwi_xy(av), 'meters',             &
2560                                     'zw(nzb_w_inner)', 424, 425, 426 )
2561
2562          ENDIF
2563
2564!
2565!--       Define the variables
2566          var_list = ';'
2567          i = 1
2568
2569          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2570
2571             IF ( INDEX( do2d(av,i), 'xy' ) /= 0 )  THEN
2572!
2573!--             If there is a star in the variable name (u* or t*), it is a
2574!--             surface variable. Define it with id_dim_zu1_xy.
2575                IF ( INDEX( do2d(av,i), '*' ) /= 0 )  THEN
2576
2577                   CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),  &
2578                                           id_dim_y_xy(av), id_dim_zu1_xy(av), &
2579                                           id_dim_time_xy(av) /), do2d(av,i),  &
2580                                           nc_precision(1), id_var_do2d(av,i), &
2581                                           TRIM( do2d_unit(av,i) ),            &
2582                                           do2d(av,i), 119, 120, 354, .TRUE. )
2583
2584                ELSE
2585
2586!
2587!--                Check for the grid
2588                   found = .FALSE.
2589                   SELECT CASE ( do2d(av,i) )
2590!
2591!--                   Most variables are defined on the zu grid
2592                      CASE ( 'e_xy', 'nc_xy', 'nr_xy', 'p_xy',                 &
2593                             'pc_xy', 'pr_xy', 'prr_xy', 'q_xy',               &
2594                             'qc_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy',           &
2595                             'ql_vp_xy', 'qr_xy', 'qv_xy',                     &
2596                             's_xy',                                           &
2597                             'theta_xy', 'thetal_xy', 'thetav_xy' )
2598
2599                         grid_x = 'x'
2600                         grid_y = 'y'
2601                         grid_z = 'zu'
2602!
2603!--                   u grid
2604                      CASE ( 'u_xy' )
2605
2606                         grid_x = 'xu'
2607                         grid_y = 'y'
2608                         grid_z = 'zu'
2609!
2610!--                   v grid
2611                      CASE ( 'v_xy' )
2612
2613                         grid_x = 'x'
2614                         grid_y = 'yv'
2615                         grid_z = 'zu'
2616!
2617!--                   w grid
2618                      CASE ( 'w_xy' )
2619
2620                         grid_x = 'x'
2621                         grid_y = 'y'
2622                         grid_z = 'zw'
2623
2624
2625                      CASE DEFAULT
2626!
2627!--                      Check for land surface quantities
2628                         IF ( land_surface )  THEN
2629                            CALL lsm_define_netcdf_grid( do2d(av,i), found,    &
2630                                                   grid_x, grid_y, grid_z )
2631                         ENDIF
2632
2633                         IF ( .NOT. found )  THEN
2634                            CALL tcm_define_netcdf_grid( do2d(av,i), found,    &
2635                                                         grid_x, grid_y,       &
2636                                                         grid_z )
2637                         ENDIF
2638
2639!
2640!--                      Check for ocean quantities
2641                         IF ( .NOT. found  .AND.  ocean_mode )  THEN
2642                            CALL ocean_define_netcdf_grid( do2d(av,i), found,  &
2643                                                           grid_x, grid_y,     &
2644                                                           grid_z )
2645                         ENDIF
2646!
2647!--                      Check for radiation quantities
2648                         IF ( .NOT. found  .AND.  radiation )  THEN
2649                            CALL radiation_define_netcdf_grid( do2d(av,i),     &
2650                                                         found, grid_x, grid_y,&
2651                                                         grid_z )
2652                         ENDIF
2653
2654!
2655!--                      Check for SALSA quantities
2656                         IF ( .NOT. found  .AND.  salsa )  THEN
2657                            CALL salsa_define_netcdf_grid( do2d(av,i), found,  &
2658                                                           grid_x, grid_y,     &
2659                                                           grid_z )
2660                         ENDIF
2661
2662!
2663!--                      Check for gust module quantities
2664                         IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
2665                            CALL gust_define_netcdf_grid( do2d(av,i), found,   &
2666                                                          grid_x, grid_y,      &
2667                                                          grid_z )
2668                         ENDIF
2669!
2670!--                      Check for biometeorology quantities
2671                         IF ( .NOT. found  .AND.  biometeorology )  THEN
2672                            CALL bio_define_netcdf_grid( do2d( av, i), found,  &
2673                                                         grid_x, grid_y,       &
2674                                                         grid_z )
2675                         ENDIF
2676!
2677!--                      Check for chemistry quantities
2678                         IF ( .NOT. found  .AND.  air_chemistry )  THEN
2679                            CALL chem_define_netcdf_grid( do2d(av,i), found,   &
2680                                                          grid_x, grid_y,      &
2681                                                          grid_z )
2682                         ENDIF
2683
2684                         IF ( .NOT. found )                                    &
2685                            CALL doq_define_netcdf_grid(                       &
2686                                                    do2d(av,i), found, grid_x, &
2687                                                    grid_y, grid_z              )
2688!
2689!--                      Check for user-defined quantities
2690                         IF ( .NOT. found  .AND.  user_module_enabled )  THEN
2691                            CALL user_define_netcdf_grid( do2d(av,i), found,   &
2692                                                          grid_x, grid_y,      &
2693                                                          grid_z )
2694                         ENDIF
2695
2696                         IF ( .NOT. found )  THEN
2697                            WRITE ( message_string, * ) 'no grid defined for', &
2698                                                ' variable ', TRIM( do2d(av,i) )
2699                            CALL message( 'define_netcdf_header', 'PA0244',    &
2700                                          0, 1, 0, 6, 0 )
2701                         ENDIF
2702
2703                   END SELECT
2704
2705!
2706!--                Select the respective dimension ids
2707                   IF ( grid_x == 'x' )  THEN
2708                      id_x = id_dim_x_xy(av)
2709                   ELSEIF ( grid_x == 'xu' )  THEN
2710                      id_x = id_dim_xu_xy(av)
2711                   ENDIF
2712
2713                   IF ( grid_y == 'y' )  THEN
2714                      id_y = id_dim_y_xy(av)
2715                   ELSEIF ( grid_y == 'yv' )  THEN
2716                      id_y = id_dim_yv_xy(av)
2717                   ENDIF
2718
2719                   IF ( grid_z == 'zu' )  THEN
2720                      id_z = id_dim_zu_xy(av)
2721                   ELSEIF ( grid_z == 'zw' )  THEN
2722                      id_z = id_dim_zw_xy(av)
2723                   ELSEIF ( grid_z == 'zs' )  THEN
2724                      id_z = id_dim_zs_xy(av)
2725                   ELSEIF ( grid_z == 'zu1' )  THEN
2726                      id_z = id_dim_zu1_xy(av)
2727                   ENDIF
2728
2729!
2730!--                Define the grid
2731                   CALL netcdf_create_var( id_set_xy(av), (/ id_x, id_y, id_z, &
2732                                           id_dim_time_xy(av) /), do2d(av,i),  &
2733                                           nc_precision(1), id_var_do2d(av,i), &
2734                                           TRIM( do2d_unit(av,i) ),            &
2735                                           do2d(av,i), 119, 120, 354, .TRUE. )
2736
2737                ENDIF
2738
2739#if defined( __netcdf4_parallel )
2740                IF ( netcdf_data_format > 4 )  THEN
2741!
2742!--                Set no fill for every variable to increase performance.
2743                   nc_stat = NF90_DEF_VAR_FILL( id_set_xy(av),     &
2744                                                id_var_do2d(av,i), &
2745                                                NF90_NOFILL, 0 )
2746                   CALL netcdf_handle_error( 'netcdf_define_header', 533 )
2747!
2748!--                Set collective io operations for parallel io
2749                   nc_stat = NF90_VAR_PAR_ACCESS( id_set_xy(av),     &
2750                                                  id_var_do2d(av,i), &
2751                                                  NF90_COLLECTIVE )
2752                   CALL netcdf_handle_error( 'netcdf_define_header', 448 )
2753                ENDIF
2754#endif
2755                var_list = TRIM( var_list) // TRIM( do2d(av,i) ) // ';'
2756
2757             ENDIF
2758
2759             i = i + 1
2760
2761          ENDDO
2762
2763!
2764!--       No arrays to output. Close the netcdf file and return.
2765          IF ( i == 1 )  RETURN
2766
2767!
2768!--       Write the list of variables as global attribute (this is used by
2769!--       restart runs and by combine_plot_fields)
2770          nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'VAR_LIST', &
2771                                  var_list )
2772          CALL netcdf_handle_error( 'netcdf_define_header', 121 )
2773
2774!
2775!--       Set general no fill, otherwise the performance drops significantly for
2776!--       parallel output.
2777          nc_stat = NF90_SET_FILL( id_set_xy(av), NF90_NOFILL, oldmode )
2778          CALL netcdf_handle_error( 'netcdf_define_header', 529 )
2779
2780!
2781!--       Leave netCDF define mode
2782          nc_stat = NF90_ENDDEF( id_set_xy(av) )
2783          CALL netcdf_handle_error( 'netcdf_define_header', 122 )
2784
2785!
2786!--       These data are only written by PE0 for parallel output to increase
2787!--       the performance.
2788          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
2789
2790!
2791!--          Write axis data: z_xy, x, y
2792             ALLOCATE( netcdf_data(1:ns) )
2793
2794!
2795!--          Write zu data
2796             DO  i = 1, ns
2797                IF( section(i,1) == -1 )  THEN
2798                   netcdf_data(i) = -1.0_wp  ! section averaged along z
2799                ELSE
2800                   netcdf_data(i) = zu( section(i,1) )
2801                ENDIF
2802             ENDDO
2803             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu_xy(av), &
2804                                     netcdf_data, start = (/ 1 /),    &
2805                                     count = (/ ns /) )
2806             CALL netcdf_handle_error( 'netcdf_define_header', 123 )
2807
2808!
2809!--          Write zw data
2810             DO  i = 1, ns
2811                IF( section(i,1) == -1 )  THEN
2812                   netcdf_data(i) = -1.0_wp  ! section averaged along z
2813                ELSE
2814                   netcdf_data(i) = zw( section(i,1) )
2815                ENDIF
2816             ENDDO
2817             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zw_xy(av), &
2818                                     netcdf_data, start = (/ 1 /),    &
2819                                     count = (/ ns /) )
2820             CALL netcdf_handle_error( 'netcdf_define_header', 124 )
2821
2822!
2823!--          Write zs data
2824             IF ( land_surface )  THEN
2825                ns_do = 0
2826                DO  i = 1, ns
2827                   IF( section(i,1) == -1 )  THEN
2828                      netcdf_data(i) = 1.0_wp  ! section averaged along z
2829                      ns_do = ns_do + 1
2830                   ELSEIF ( section(i,1) < nzs )  THEN
2831                      netcdf_data(i) = - zs( section(i,1) )
2832                      ns_do = ns_do + 1
2833                   ENDIF
2834                ENDDO
2835
2836                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zs_xy(av), &
2837                                        netcdf_data(1:ns_do), start = (/ 1 /),    &
2838                                        count = (/ ns_do /) )
2839                CALL netcdf_handle_error( 'netcdf_define_header', 124 )
2840
2841             ENDIF
2842
2843!
2844!--          Write gridpoint number data
2845             netcdf_data(1:ns) = section(1:ns,1)
2846             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_ind_z_xy(av), &
2847                                     netcdf_data, start = (/ 1 /),       &
2848                                     count = (/ ns /) )
2849             CALL netcdf_handle_error( 'netcdf_define_header', 125 )
2850
2851             DEALLOCATE( netcdf_data )
2852
2853!
2854!--          Write the cross section height u*, t*
2855             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu1_xy(av), &
2856                                     (/ zu(nzb+1) /), start = (/ 1 /), &
2857                                     count = (/ 1 /) )
2858             CALL netcdf_handle_error( 'netcdf_define_header', 126 )
2859
2860!
2861!--          Write data for x (shifted by +dx/2) and xu axis
2862             ALLOCATE( netcdf_data(0:nx) )
2863
2864             DO  i = 0, nx
2865                netcdf_data(i) = ( i + 0.5_wp ) * dx
2866             ENDDO
2867
2868             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_x_xy(av), &
2869                                     netcdf_data, start = (/ 1 /),   &
2870                                     count = (/ nx+1 /) )
2871             CALL netcdf_handle_error( 'netcdf_define_header', 127 )
2872
2873             DO  i = 0, nx
2874                netcdf_data(i) = i * dx
2875             ENDDO
2876
2877             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_xu_xy(av), &
2878                                     netcdf_data, start = (/ 1 /),    &
2879                                     count = (/ nx+1 /) )
2880             CALL netcdf_handle_error( 'netcdf_define_header', 367 )
2881
2882             DEALLOCATE( netcdf_data )
2883
2884!
2885!--          Write data for y (shifted by +dy/2) and yv axis
2886             ALLOCATE( netcdf_data(0:ny+1) )
2887
2888             DO  i = 0, ny
2889                netcdf_data(i) = ( i + 0.5_wp ) * dy
2890             ENDDO
2891
2892             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_y_xy(av), &
2893                                     netcdf_data, start = (/ 1 /),   &
2894                                     count = (/ ny+1 /))
2895             CALL netcdf_handle_error( 'netcdf_define_header', 128 )
2896
2897             DO  i = 0, ny
2898                netcdf_data(i) = i * dy
2899             ENDDO
2900
2901             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_yv_xy(av), &
2902                                     netcdf_data, start = (/ 1 /),    &
2903                                     count = (/ ny+1 /))
2904             CALL netcdf_handle_error( 'netcdf_define_header', 368 )
2905
2906             DEALLOCATE( netcdf_data )
2907!
2908!--          Write UTM coordinates
2909             IF ( init_model%rotation_angle == 0.0_wp )  THEN
2910!
2911!--             1D in case of no rotation
2912                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
2913!
2914!--             x coordinates
2915                ALLOCATE( netcdf_data(0:nx) )
2916                DO  k = 0, 2
2917!
2918!--                Scalar grid points
2919                   IF ( k == 0 )  THEN
2920                      shift_x = 0.5
2921!
2922!--                u grid points
2923                   ELSEIF ( k == 1 )  THEN
2924                      shift_x = 0.0
2925!
2926!--                v grid points
2927                   ELSEIF ( k == 2 )  THEN
2928                      shift_x = 0.5
2929                   ENDIF
2930
2931                   DO  i = 0, nx
2932                     netcdf_data(i) = init_model%origin_x            &
2933                                    + cos_ra * ( i + shift_x ) * dx
2934                   ENDDO
2935
2936                   nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_eutm_xy(k,av),&
2937                                           netcdf_data, start = (/ 1 /),   &
2938                                           count = (/ nx+1 /) )
2939                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
2940
2941                ENDDO
2942                DEALLOCATE( netcdf_data )
2943!
2944!--             y coordinates
2945                ALLOCATE( netcdf_data(0:ny) )
2946                DO  k = 0, 2
2947!
2948!--                Scalar grid points
2949                   IF ( k == 0 )  THEN
2950                      shift_y = 0.5
2951!
2952!--                u grid points
2953                   ELSEIF ( k == 1 )  THEN
2954                      shift_y = 0.5
2955!
2956!--                v grid points
2957                   ELSEIF ( k == 2 )  THEN
2958                      shift_y = 0.0
2959                   ENDIF
2960
2961                   DO  j = 0, ny
2962                      netcdf_data(j) = init_model%origin_y            &
2963                                     + cos_ra * ( j + shift_y ) * dy
2964                   ENDDO
2965
2966                   nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_nutm_xy(k,av),&
2967                                           netcdf_data, start = (/ 1 /),   &
2968                                           count = (/ ny+1 /) )
2969                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
2970
2971                ENDDO
2972                DEALLOCATE( netcdf_data )
2973
2974             ELSE
2975!
2976!--             2D in case of rotation
2977                ALLOCATE( netcdf_data_2d(0:nx,0:ny) )
2978                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
2979                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
2980
2981                DO  k = 0, 2
2982!
2983!--               Scalar grid points
2984                  IF ( k == 0 )  THEN
2985                     shift_x = 0.5 ; shift_y = 0.5
2986!
2987!--               u grid points
2988                  ELSEIF ( k == 1 )  THEN
2989                     shift_x = 0.0 ; shift_y = 0.5
2990!
2991!--               v grid points
2992                  ELSEIF ( k == 2 )  THEN
2993                     shift_x = 0.5 ; shift_y = 0.0
2994                  ENDIF
2995
2996                  DO  j = 0, ny
2997                     DO  i = 0, nx
2998                        netcdf_data_2d(i,j) = init_model%origin_x            &
2999                                            + cos_ra * ( i + shift_x ) * dx  &
3000                                            + sin_ra * ( j + shift_y ) * dy
3001                     ENDDO
3002                  ENDDO
3003
3004                  nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_eutm_xy(k,av),  &
3005                                          netcdf_data_2d, start = (/ 1, 1 /),   &
3006                                          count = (/ nx+1, ny+1 /) )
3007                  CALL netcdf_handle_error( 'netcdf_define_header', 555 )
3008
3009                  DO  j = 0, ny
3010                     DO  i = 0, nx
3011                        netcdf_data_2d(i,j) = init_model%origin_y            &
3012                                            - sin_ra * ( i + shift_x ) * dx  &
3013                                            + cos_ra * ( j + shift_y ) * dy
3014                     ENDDO
3015                  ENDDO
3016
3017                  nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_nutm_xy(k,av),  &
3018                                          netcdf_data_2d, start = (/ 1, 1 /),   &
3019                                          count = (/ nx+1, ny+1 /) )
3020                  CALL netcdf_handle_error( 'netcdf_define_header', 556 )
3021
3022                ENDDO
3023                DEALLOCATE( netcdf_data_2d )
3024             ENDIF
3025
3026          ENDIF
3027!
3028!--       Write lon and lat data. Only for parallel output.
3029          IF ( netcdf_data_format > 4 )  THEN
3030
3031             ALLOCATE( lat(nxl:nxr,nys:nyn) )
3032             ALLOCATE( lon(nxl:nxr,nys:nyn) )
3033             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
3034             sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
3035
3036             DO  k = 0, 2
3037!
3038!--             Scalar grid points
3039                IF ( k == 0 )  THEN
3040                   shift_x = 0.5 ; shift_y = 0.5
3041!
3042!--             u grid points
3043                ELSEIF ( k == 1 )  THEN
3044                   shift_x = 0.0 ; shift_y = 0.5
3045!
3046!--             v grid points
3047                ELSEIF ( k == 2 )  THEN
3048                   shift_x = 0.5 ; shift_y = 0.0
3049                ENDIF
3050
3051                DO  j = nys, nyn
3052                   DO  i = nxl, nxr
3053                      eutm = init_model%origin_x            &
3054                           + cos_ra * ( i + shift_x ) * dx  &
3055                           + sin_ra * ( j + shift_y ) * dy
3056                      nutm = init_model%origin_y            &
3057                           - sin_ra * ( i + shift_x ) * dx  &
3058                           + cos_ra * ( j + shift_y ) * dy
3059
3060                      CALL  convert_utm_to_geographic( crs_list,          &
3061                                                       eutm, nutm,        &
3062                                                       lon(i,j), lat(i,j) )
3063                   ENDDO
3064                ENDDO
3065
3066                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_lon_xy(k,av), &
3067                                     lon, start = (/ nxl+1, nys+1 /),       &
3068                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
3069                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
3070
3071                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_lat_xy(k,av), &
3072                                     lat, start = (/ nxl+1, nys+1 /),       &
3073                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
3074                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
3075             ENDDO
3076
3077             DEALLOCATE( lat )
3078             DEALLOCATE( lon )
3079
3080          ENDIF
3081!
3082!--       In case of non-flat topography write height information. Only for
3083!--       parallel netcdf output.
3084          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
3085               netcdf_data_format > 4  )  THEN
3086
3087!             IF ( nxr == nx  .AND.  nyn /= ny )  THEN
3088!                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av),     &
3089!                                        zu_s_inner(nxl:nxr+1,nys:nyn),         &
3090!                                        start = (/ nxl+1, nys+1 /),            &
3091!                                        count = (/ nxr-nxl+2, nyn-nys+1 /) )
3092!             ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
3093!                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av),     &
3094!                                        zu_s_inner(nxl:nxr,nys:nyn+1),         &
3095!                                        start = (/ nxl+1, nys+1 /),            &
3096!                                        count = (/ nxr-nxl+1, nyn-nys+2 /) )
3097!             ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
3098!                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av),     &
3099!                                        zu_s_inner(nxl:nxr+1,nys:nyn+1),       &
3100!                                        start = (/ nxl+1, nys+1 /),            &
3101!                                        count = (/ nxr-nxl+2, nyn-nys+2 /) )
3102!             ELSE
3103                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av),     &
3104                                        zu_s_inner(nxl:nxr,nys:nyn),           &
3105                                        start = (/ nxl+1, nys+1 /),            &
3106                                        count = (/ nxr-nxl+1, nyn-nys+1 /) )
3107!             ENDIF
3108             CALL netcdf_handle_error( 'netcdf_define_header', 427 )
3109
3110!             IF ( nxr == nx  .AND.  nyn /= ny )  THEN
3111!                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av),     &
3112!                                        zw_w_inner(nxl:nxr+1,nys:nyn),         &
3113!                                        start = (/ nxl+1, nys+1 /),            &
3114!                                        count = (/ nxr-nxl+2, nyn-nys+1 /) )
3115!             ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
3116!                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av),     &
3117!                                        zw_w_inner(nxl:nxr,nys:nyn+1),         &
3118!                                        start = (/ nxl+1, nys+1 /),            &
3119!                                        count = (/ nxr-nxl+1, nyn-nys+2 /) )
3120!             ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
3121!                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av),     &
3122!                                        zw_w_inner(nxl:nxr+1,nys:nyn+1),       &
3123!                                        start = (/ nxl+1, nys+1 /),            &
3124!                                        count = (/ nxr-nxl+2, nyn-nys+2 /) )
3125!             ELSE
3126                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av),     &
3127                                        zw_w_inner(nxl:nxr,nys:nyn),           &
3128                                        start = (/ nxl+1, nys+1 /),            &
3129                                        count = (/ nxr-nxl+1, nyn-nys+1 /) )
3130!             ENDIF
3131             CALL netcdf_handle_error( 'netcdf_define_header', 428 )
3132
3133          ENDIF
3134
3135       CASE ( 'xy_ext' )
3136
3137!
3138!--       Get the list of variables and compare with the actual run.
3139!--       First var_list_old has to be reset, since GET_ATT does not assign
3140!--       trailing blanks.
3141          var_list_old = ' '
3142          nc_stat = NF90_GET_ATT( id_set_xy(av), NF90_GLOBAL, 'VAR_LIST', &
3143                                  var_list_old )
3144          CALL netcdf_handle_error( 'netcdf_define_header', 129 )
3145
3146          var_list = ';'
3147          i = 1
3148          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
3149             IF ( INDEX( do2d(av,i), 'xy' ) /= 0 )  THEN
3150                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
3151             ENDIF
3152             i = i + 1
3153          ENDDO
3154
3155          IF ( av == 0 )  THEN
3156             var = '(xy)'
3157          ELSE
3158             var = '(xy_av)'
3159          ENDIF
3160
3161          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
3162             message_string = 'netCDF file for cross-sections ' //           &
3163                              TRIM( var ) // ' from previous run found,' //  &
3164                              '&but this file cannot be extended due to' //  &
3165                              ' variable mismatch.' //                       &
3166                              '&New file is created instead.'
3167             CALL message( 'define_netcdf_header', 'PA0249', 0, 1, 0, 6, 0 )
3168             extend = .FALSE.
3169             RETURN
3170          ENDIF
3171
3172!
3173!--       Calculate the number of current sections
3174          ns = 1
3175          DO WHILE ( section(ns,1) /= -9999  .AND.  ns <= 100 )
3176             ns = ns + 1
3177          ENDDO
3178          ns = ns - 1
3179
3180!
3181!--       Get and compare the number of horizontal cross sections
3182          nc_stat = NF90_INQ_VARID( id_set_xy(av), 'zu_xy', id_var_zu_xy(av) )
3183          CALL netcdf_handle_error( 'netcdf_define_header', 130 )
3184
3185          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_zu_xy(av), &
3186                                           dimids = id_dim_zu_xy_old )
3187          CALL netcdf_handle_error( 'netcdf_define_header', 131 )
3188          id_dim_zu_xy(av) = id_dim_zu_xy_old(1)
3189
3190          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_zu_xy(av), &
3191                                            len = ns_old )
3192          CALL netcdf_handle_error( 'netcdf_define_header', 132 )
3193
3194          IF ( ns /= ns_old )  THEN
3195             message_string = 'netCDF file for cross-sections ' //          &
3196                              TRIM( var ) // ' from previous run found,' // &
3197                              '&but this file cannot be extended due to' // &
3198                              ' mismatch in number of' //                   &
3199                              ' cross sections.' //                         &
3200                              '&New file is created instead.'
3201             CALL message( 'define_netcdf_header', 'PA0250', 0, 1, 0, 6, 0 )
3202             extend = .FALSE.
3203             RETURN
3204          ENDIF
3205
3206!
3207!--       Get and compare the heights of the cross sections
3208          ALLOCATE( netcdf_data(1:ns_old) )
3209
3210          nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_zu_xy(av), netcdf_data )
3211          CALL netcdf_handle_error( 'netcdf_define_header', 133 )
3212
3213          DO  i = 1, ns
3214             IF ( section(i,1) /= -1 )  THEN
3215                IF ( zu(section(i,1)) /= netcdf_data(i) )  THEN
3216                   message_string = 'netCDF file for cross-sections ' //       &
3217                               TRIM( var ) // ' from previous run found,' //   &
3218                               ' but this file cannot be extended' //          &
3219                               ' due to mismatch in cross' //                  &
3220                               ' section levels.' //                           &
3221                               ' New file is created instead.'
3222                   CALL message( 'define_netcdf_header', 'PA0251',             &
3223                                                                 0, 1, 0, 6, 0 )
3224                   extend = .FALSE.
3225                   RETURN
3226                ENDIF
3227             ELSE
3228                IF ( -1.0_wp /= netcdf_data(i) )  THEN
3229                   message_string = 'netCDF file for cross-sections ' //       &
3230                               TRIM( var ) // ' from previous run found,' //   &
3231                               ' but this file cannot be extended' //          &
3232                               ' due to mismatch in cross' //                  &
3233                               ' section levels.' //                           &
3234                               ' New file is created instead.'
3235                   CALL message( 'define_netcdf_header', 'PA0251',             &
3236                                                                 0, 1, 0, 6, 0 )
3237                   extend = .FALSE.
3238                   RETURN
3239                ENDIF
3240             ENDIF
3241          ENDDO
3242
3243          DEALLOCATE( netcdf_data )
3244
3245!
3246!--       Get the id of the time coordinate (unlimited coordinate) and its
3247!--       last index on the file. The next time level is do2d..count+1.
3248!--       The current time must be larger than the last output time
3249!--       on the file.
3250          nc_stat = NF90_INQ_VARID( id_set_xy(av), 'time', id_var_time_xy(av) )
3251          CALL netcdf_handle_error( 'netcdf_define_header', 134 )
3252
3253          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_time_xy(av), &
3254                                           dimids = id_dim_time_old )
3255          CALL netcdf_handle_error( 'netcdf_define_header', 135 )
3256          id_dim_time_xy(av) = id_dim_time_old(1)
3257
3258          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_time_xy(av), &
3259                                            len = ntime_count )
3260          CALL netcdf_handle_error( 'netcdf_define_header', 136 )
3261
3262!
3263!--       For non-parallel output use the last output time level of the netcdf
3264!--       file because the time dimension is unlimited. In case of parallel
3265!--       output the variable ntime_count could get the value of 9*10E36 because
3266!--       the time dimension is limited.
3267          IF ( netcdf_data_format < 5 ) do2d_xy_time_count(av) = ntime_count
3268
3269          nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_time_xy(av),           &
3270                                  last_time_coordinate,                        &
3271                                  start = (/ do2d_xy_time_count(av) /),        &
3272                                  count = (/ 1 /) )
3273          CALL netcdf_handle_error( 'netcdf_define_header', 137 )
3274
3275          IF ( last_time_coordinate(1) >= simulated_time )  THEN
3276             message_string = 'netCDF file for cross sections ' //             &
3277                              TRIM( var ) // ' from previous run found,' //    &
3278                              '&but this file cannot be extended becaus' //    &
3279                              'e the current output time' //                   &
3280                              '&is less or equal than the last output t' //    &
3281                              'ime on this file.' //                           &
3282                              '&New file is created instead.'
3283             CALL message( 'define_netcdf_header', 'PA0252', 0, 1, 0, 6, 0 )
3284             do2d_xy_time_count(av) = 0
3285             extend = .FALSE.
3286             RETURN
3287          ENDIF
3288
3289          IF ( netcdf_data_format > 4 )  THEN
3290!
3291!--          Check if the needed number of output time levels is increased
3292!--          compared to the number of time levels in the existing file.
3293             IF ( ntdim_2d_xy(av) > ntime_count )  THEN
3294                message_string = 'netCDF file for cross sections ' //          &
3295                                 TRIM( var ) // ' from previous run found,' // &
3296                                 '&but this file cannot be extended becaus' // &
3297                                 'e the number of output time levels has b' // &
3298                                 'een increased compared to the previous s' // &
3299                                 'imulation.' //                               &
3300                                 '&New file is created instead.'
3301                CALL message( 'define_netcdf_header', 'PA0389', 0, 1, 0, 6, 0 )
3302                do2d_xy_time_count(av) = 0
3303                extend = .FALSE.
3304!
3305!--             Recalculate the needed time levels for the new file.
3306                IF ( av == 0 )  THEN
3307                   ntdim_2d_xy(0) = CEILING(                            &
3308                           ( end_time - MAX( skip_time_do2d_xy,         &
3309                                             simulated_time_at_begin )  &
3310                           ) / dt_do2d_xy )
3311                   IF ( do2d_at_begin )  ntdim_2d_xy(0) = ntdim_2d_xy(0) + 1
3312                ELSE
3313                   ntdim_2d_xy(1) = CEILING(                            &
3314                           ( end_time - MAX( skip_time_data_output_av,  &
3315                                             simulated_time_at_begin )  &
3316                           ) / dt_data_output_av )
3317                ENDIF
3318                RETURN
3319             ENDIF
3320          ENDIF
3321
3322!
3323!--       Dataset seems to be extendable.
3324!--       Now get the variable ids.
3325          i = 1
3326          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
3327             IF ( INDEX( do2d(av,i), 'xy' ) /= 0 )  THEN
3328                nc_stat = NF90_INQ_VARID( id_set_xy(av), do2d(av,i), &
3329                                          id_var_do2d(av,i) )
3330                CALL netcdf_handle_error( 'netcdf_define_header', 138 )
3331#if defined( __netcdf4_parallel )
3332!
3333!--             Set collective io operations for parallel io
3334                IF ( netcdf_data_format > 4 )  THEN
3335                   nc_stat = NF90_VAR_PAR_ACCESS( id_set_xy(av),     &
3336                                                  id_var_do2d(av,i), &
3337                                                  NF90_COLLECTIVE )
3338                   CALL netcdf_handle_error( 'netcdf_define_header', 454 )
3339                ENDIF
3340#endif
3341             ENDIF
3342             i = i + 1
3343          ENDDO
3344
3345!
3346!--       Update the title attribute on file
3347!--       In order to avoid 'data mode' errors if updated attributes are larger
3348!--       than their original size, NF90_PUT_ATT is called in 'define mode'
3349!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
3350!--       performance loss due to data copying; an alternative strategy would be
3351!--       to ensure equal attribute size in a job chain. Maybe revise later.
3352          IF ( av == 0 )  THEN
3353             time_average_text = ' '
3354          ELSE
3355             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
3356                                                            averaging_interval
3357          ENDIF
3358          nc_stat = NF90_REDEF( id_set_xy(av) )
3359          CALL netcdf_handle_error( 'netcdf_define_header', 431 )
3360          nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'title',         &
3361                                  TRIM( run_description_header ) //            &
3362                                  TRIM( time_average_text ) )
3363          CALL netcdf_handle_error( 'netcdf_define_header', 139 )
3364          nc_stat = NF90_ENDDEF( id_set_xy(av) )
3365          CALL netcdf_handle_error( 'netcdf_define_header', 432 )
3366          message_string = 'netCDF file for cross-sections ' //                &
3367                            TRIM( var ) // ' from previous run found.' //      &
3368                           '&This file will be extended.'
3369          CALL message( 'define_netcdf_header', 'PA0253', 0, 0, 0, 6, 0 )
3370
3371
3372       CASE ( 'xz_new' )
3373
3374!
3375!--       Define some global attributes of the dataset
3376          IF ( av == 0 )  THEN
3377             CALL netcdf_create_global_atts( id_set_xz(av), 'xz', TRIM( run_description_header ), 140 )
3378             time_average_text = ' '
3379          ELSE
3380             CALL netcdf_create_global_atts( id_set_xz(av), 'xz_av', TRIM( run_description_header ), 140 )
3381             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
3382             nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'time_avg',   &
3383                                     TRIM( time_average_text ) )
3384             CALL netcdf_handle_error( 'netcdf_define_header', 141 )
3385          ENDIF
3386
3387!
3388!--       Define time coordinate for xz sections.
3389!--       For parallel output the time dimensions has to be limited, otherwise
3390!--       the performance drops significantly.
3391          IF ( netcdf_data_format < 5 )  THEN
3392             CALL netcdf_create_dim( id_set_xz(av), 'time', NF90_UNLIMITED,    &
3393                                     id_dim_time_xz(av), 142 )
3394          ELSE
3395             CALL netcdf_create_dim( id_set_xz(av), 'time', ntdim_2d_xz(av),   &
3396                                     id_dim_time_xz(av), 525 )
3397          ENDIF
3398
3399          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_time_xz(av) /),     &
3400                                  'time', NF90_DOUBLE, id_var_time_xz(av),     &
3401                                  'seconds', 'time', 143, 144, 000 )
3402          CALL netcdf_create_att( id_set_xz(av), id_var_time_xz(av), 'standard_name', 'time', 000)
3403          CALL netcdf_create_att( id_set_xz(av), id_var_time_xz(av), 'axis', 'T', 000)
3404!
3405!--       Define the spatial dimensions and coordinates for xz-sections.
3406!--       First, determine the number of vertical sections to be written.
3407          IF ( section(1,2) == -9999 )  THEN
3408             RETURN
3409          ELSE
3410             ns = 1
3411             DO WHILE ( section(ns,2) /= -9999  .AND.  ns <= 100 )
3412                ns = ns + 1
3413             ENDDO
3414             ns = ns - 1
3415          ENDIF
3416
3417!
3418!--       Define y-axis (for scalar position)
3419          CALL netcdf_create_dim( id_set_xz(av), 'y_xz', ns, id_dim_y_xz(av),  &
3420                                  145 )
3421          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_y_xz(av) /),        &
3422                                  'y_xz', NF90_DOUBLE, id_var_y_xz(av),        &
3423                                  'meters', '', 146, 147, 000 )
3424!
3425!--       Define y-axis (for v position)
3426          CALL netcdf_create_dim( id_set_xz(av), 'yv_xz', ns,                  &
3427                                  id_dim_yv_xz(av), 369 )
3428          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_yv_xz(av) /),       &
3429                                  'yv_xz', NF90_DOUBLE, id_var_yv_xz(av),      &
3430                                  'meters', '', 370, 371, 000 )
3431!
3432!--       Define a variable to store the layer indices of the vertical cross
3433!--       sections
3434          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_y_xz(av) /),        &
3435                                  'ind_y_xz', NF90_DOUBLE,                     &
3436                                  id_var_ind_y_xz(av), 'gridpoints', '', 148,  &
3437                                  149, 000 )
3438!
3439!--       Define x-axis (for scalar position)
3440          CALL netcdf_create_dim( id_set_xz(av), 'x', nx+1, id_dim_x_xz(av),   &
3441                                  150 )
3442          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_x_xz(av) /), 'x',   &
3443                                  NF90_DOUBLE, id_var_x_xz(av), 'meters', '',  &
3444                                  151, 152, 000 )
3445!
3446!--       Define x-axis (for u position)
3447          CALL netcdf_create_dim( id_set_xz(av), 'xu', nx+1, id_dim_xu_xz(av), &
3448                                  372 )
3449          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_xu_xz(av) /), 'xu', &
3450                                  NF90_DOUBLE, id_var_xu_xz(av), 'meters', '', &
3451                                  373, 374, 000 )
3452!
3453!--       Define the three z-axes (zu, zw, and zs)
3454          CALL netcdf_create_dim( id_set_xz(av), 'zu', nz+2, id_dim_zu_xz(av), &
3455                                  153 )
3456          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zu_xz(av) /), 'zu', &
3457                                  NF90_DOUBLE, id_var_zu_xz(av), 'meters', '', &
3458                                  154, 155, 000 )
3459          CALL netcdf_create_dim( id_set_xz(av), 'zw', nz+2, id_dim_zw_xz(av), &
3460                                  156 )
3461          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zw_xz(av) /), 'zw', &
3462                                  NF90_DOUBLE, id_var_zw_xz(av), 'meters', '', &
3463                                  157, 158, 000 )
3464!
3465!--       Define UTM and geographic coordinates
3466          CALL define_geo_coordinates( id_set_xz(av),         &
3467                  (/ id_dim_x_xz(av), id_dim_xu_xz(av) /),    &
3468                  (/ id_dim_y_xz(av), id_dim_yv_xz(av) /),    &
3469                  id_var_eutm_xz(:,av), id_var_nutm_xz(:,av), &
3470                  id_var_lat_xz(:,av), id_var_lon_xz(:,av)    )
3471!
3472!--       Define coordinate-reference system
3473          CALL netcdf_create_crs( id_set_xz(av), 000 )
3474
3475          IF ( land_surface )  THEN
3476
3477             CALL netcdf_create_dim( id_set_xz(av), 'zs', nzs,                 &
3478                                     id_dim_zs_xz(av), 542 )
3479             CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zs_xz(av) /),    &
3480                                     'zs', NF90_DOUBLE, id_var_zs_xz(av),      &
3481                                     'meters', '', 543, 544, 000 )
3482
3483          ENDIF
3484
3485!
3486!--       Define the variables
3487          var_list = ';'
3488          i = 1
3489
3490          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
3491
3492             IF ( INDEX( do2d(av,i), 'xz' ) /= 0 )  THEN
3493
3494!
3495!--             Check for the grid
3496                found = .FALSE.
3497                SELECT CASE ( do2d(av,i) )
3498!
3499!--                Most variables are defined on the zu grid
3500                   CASE ( 'e_xz', 'nc_xz', 'nr_xz', 'p_xz', 'pc_xz',           &
3501                          'pr_xz', 'prr_xz', 'q_xz', 'qc_xz',                  &
3502                          'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz',  &
3503                          'qv_xz', 's_xz',                                     &
3504                          'theta_xz', 'thetal_xz', 'thetav_xz' )
3505
3506                      grid_x = 'x'
3507                      grid_y = 'y'
3508                      grid_z = 'zu'
3509!
3510!--                u grid
3511                   CASE ( 'u_xz' )
3512
3513                      grid_x = 'xu'
3514                      grid_y = 'y'
3515                      grid_z = 'zu'
3516!
3517!--                v grid
3518                   CASE ( 'v_xz' )
3519
3520                      grid_x = 'x'
3521                      grid_y = 'yv'
3522                      grid_z = 'zu'
3523!
3524!--                w grid
3525                   CASE ( 'w_xz' )
3526
3527                      grid_x = 'x'
3528                      grid_y = 'y'
3529                      grid_z = 'zw'
3530
3531                   CASE DEFAULT
3532
3533!
3534!--                   Check for land surface quantities
3535                      IF ( land_surface )  THEN
3536                         CALL lsm_define_netcdf_grid( do2d(av,i), found,       &
3537                                                      grid_x, grid_y, grid_z )
3538                      ENDIF
3539
3540                      IF ( .NOT. found )  THEN
3541                         CALL tcm_define_netcdf_grid( do2d(av,i), found,       &
3542                                                      grid_x, grid_y, grid_z )
3543                      ENDIF
3544
3545!
3546!--                   Check for ocean quantities
3547                      IF ( .NOT. found  .AND.  ocean_mode )  THEN
3548                         CALL ocean_define_netcdf_grid( do2d(av,i), found,  &
3549                                                        grid_x, grid_y, grid_z )
3550                      ENDIF
3551!
3552!--                   Check for radiation quantities
3553                      IF ( .NOT. found  .AND.  radiation )  THEN
3554                         CALL radiation_define_netcdf_grid( do2d(av,i), found, &
3555                                                            grid_x, grid_y,    &
3556                                                            grid_z )
3557                      ENDIF
3558!
3559!--                   Check for SALSA quantities
3560                      IF ( .NOT. found  .AND.  salsa )  THEN
3561                         CALL salsa_define_netcdf_grid( do2d(av,i), found,     &
3562                                                        grid_x, grid_y, grid_z )
3563                      ENDIF
3564
3565!
3566!--                   Check for gust module quantities
3567                      IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
3568                         CALL gust_define_netcdf_grid( do2d(av,i), found,      &
3569                                                       grid_x, grid_y, grid_z )
3570                      ENDIF
3571
3572!
3573!--                   Check for chemistry quantities
3574                      IF ( .NOT. found  .AND.  air_chemistry )  THEN
3575                         CALL chem_define_netcdf_grid( do2d(av,i), found,      &
3576                                                       grid_x, grid_y,         &
3577                                                       grid_z )
3578                      ENDIF
3579
3580                      IF ( .NOT. found )                                       &
3581                         CALL doq_define_netcdf_grid( do2d(av,i), found,       &
3582                                                      grid_x, grid_y, grid_z )
3583
3584!
3585!--                   Check for user-defined quantities
3586                      IF ( .NOT. found  .AND.  user_module_enabled )  THEN
3587                         CALL user_define_netcdf_grid( do2d(av,i), found,      &
3588                                                       grid_x, grid_y, grid_z )
3589                      ENDIF
3590
3591                      IF ( .NOT. found )  THEN
3592                         WRITE ( message_string, * ) 'no grid defined for',    &
3593                                                ' variable ', TRIM( do2d(av,i) )
3594                         CALL message( 'define_netcdf_header', 'PA0244',       &
3595                                       0, 1, 0, 6, 0 )
3596                      ENDIF
3597
3598                END SELECT
3599
3600!
3601!--             Select the respective dimension ids
3602                IF ( grid_x == 'x' )  THEN
3603                   id_x = id_dim_x_xz(av)
3604                ELSEIF ( grid_x == 'xu' )  THEN
3605                   id_x = id_dim_xu_xz(av)
3606                ENDIF
3607
3608                IF ( grid_y == 'y' )  THEN
3609                   id_y = id_dim_y_xz(av)
3610                ELSEIF ( grid_y == 'yv' )  THEN
3611                   id_y = id_dim_yv_xz(av)
3612                ENDIF
3613
3614                IF ( grid_z == 'zu' )  THEN
3615                   id_z = id_dim_zu_xz(av)
3616                ELSEIF ( grid_z == 'zw' )  THEN
3617                   id_z = id_dim_zw_xz(av)
3618                ELSEIF ( grid_z == 'zs' )  THEN
3619                   id_z = id_dim_zs_xz(av)
3620                ENDIF
3621
3622!
3623!--             Define the grid
3624                CALL netcdf_create_var( id_set_xz(av), (/ id_x, id_y, id_z,    &
3625                                        id_dim_time_xz(av) /), do2d(av,i),     &
3626                                        nc_precision(2), id_var_do2d(av,i),    &
3627                                        TRIM( do2d_unit(av,i) ), do2d(av,i),   &
3628                                        159, 160, 355, .TRUE. )
3629
3630#if defined( __netcdf4_parallel )
3631
3632                IF ( netcdf_data_format > 4 )  THEN
3633!
3634!--                Set no fill for every variable to increase performance.
3635                   nc_stat = NF90_DEF_VAR_FILL( id_set_xz(av),     &
3636                                                id_var_do2d(av,i), &
3637                                                NF90_NOFILL, 0 )
3638                   CALL netcdf_handle_error( 'netcdf_define_header', 534 )
3639!
3640!--                Set independent io operations for parallel io. Collective io
3641!--                is only allowed in case of a 1d-decomposition along x,
3642!--                because otherwise, not all PEs have output data.
3643                   IF ( npey == 1 )  THEN
3644                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
3645                                                     id_var_do2d(av,i), &
3646                                                     NF90_COLLECTIVE )
3647                   ELSE
3648!
3649!--                   Test simulations showed that the output of cross sections
3650!--                   by all PEs in data_output_2d using NF90_COLLECTIVE is
3651!--                   faster than the output by the first row of PEs in
3652!--                   x-direction using NF90_INDEPENDENT.
3653                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),    &
3654                                                    id_var_do2d(av,i), &
3655                                                    NF90_COLLECTIVE )
3656!                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
3657!                                                     id_var_do2d(av,i), &
3658!                                                     NF90_INDEPENDENT )
3659                   ENDIF
3660                   CALL netcdf_handle_error( 'netcdf_define_header', 449 )
3661                ENDIF
3662#endif
3663                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
3664
3665             ENDIF
3666
3667             i = i + 1
3668
3669          ENDDO
3670
3671!
3672!--       No arrays to output. Close the netcdf file and return.
3673          IF ( i == 1 )  RETURN
3674
3675!
3676!--       Write the list of variables as global attribute (this is used by
3677!--       restart runs and by combine_plot_fields)
3678          nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'VAR_LIST', &
3679                                  var_list )
3680          CALL netcdf_handle_error( 'netcdf_define_header', 161 )
3681
3682!
3683!--       Set general no fill, otherwise the performance drops significantly for
3684!--       parallel output.
3685          nc_stat = NF90_SET_FILL( id_set_xz(av), NF90_NOFILL, oldmode )
3686          CALL netcdf_handle_error( 'netcdf_define_header', 530 )
3687
3688!
3689!--       Leave netCDF define mode
3690          nc_stat = NF90_ENDDEF( id_set_xz(av) )
3691          CALL netcdf_handle_error( 'netcdf_define_header', 162 )
3692
3693!
3694!--       These data are only written by PE0 for parallel output to increase
3695!--       the performance.
3696          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
3697
3698!
3699!--          Write axis data: y_xz, x, zu, zw
3700             ALLOCATE( netcdf_data(1:ns) )
3701
3702!
3703!--          Write y_xz data (shifted by +dy/2)
3704             DO  i = 1, ns
3705                IF( section(i,2) == -1 )  THEN
3706                   netcdf_data(i) = -1.0_wp  ! section averaged along y
3707                ELSE
3708                   netcdf_data(i) = ( section(i,2) + 0.5_wp ) * dy
3709                ENDIF
3710             ENDDO
3711             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_y_xz(av), &
3712                                     netcdf_data, start = (/ 1 /),   &
3713                                     count = (/ ns /) )
3714             CALL netcdf_handle_error( 'netcdf_define_header', 163 )
3715
3716!
3717!--          Write yv_xz data
3718             DO  i = 1, ns
3719                IF( section(i,2) == -1 )  THEN
3720                   netcdf_data(i) = -1.0_wp  ! section averaged along y
3721                ELSE
3722                   netcdf_data(i) = section(i,2) * dy
3723                ENDIF
3724             ENDDO
3725             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_yv_xz(av), &
3726                                     netcdf_data, start = (/ 1 /),    &
3727                                     count = (/ ns /) )
3728             CALL netcdf_handle_error( 'netcdf_define_header', 375 )
3729
3730!
3731!--          Write gridpoint number data
3732             netcdf_data(1:ns) = section(1:ns,2)
3733             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_ind_y_xz(av), &
3734                                     netcdf_data, start = (/ 1 /),       &
3735                                     count = (/ ns /) )
3736             CALL netcdf_handle_error( 'netcdf_define_header', 164 )
3737
3738
3739             DEALLOCATE( netcdf_data )
3740
3741!
3742!--          Write data for x (shifted by +dx/2) and xu axis
3743             ALLOCATE( netcdf_data(0:nx) )
3744
3745             DO  i = 0, nx
3746                netcdf_data(i) = ( i + 0.5_wp ) * dx
3747             ENDDO
3748
3749             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_x_xz(av), &
3750                                     netcdf_data, start = (/ 1 /),   &
3751                                     count = (/ nx+1 /) )
3752             CALL netcdf_handle_error( 'netcdf_define_header', 165 )
3753
3754             DO  i = 0, nx
3755                netcdf_data(i) = i * dx
3756             ENDDO
3757
3758             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_xu_xz(av), &
3759                                     netcdf_data, start = (/ 1 /),    &
3760                                     count = (/ nx+1 /) )
3761             CALL netcdf_handle_error( 'netcdf_define_header', 377 )
3762
3763             DEALLOCATE( netcdf_data )
3764
3765!
3766!--          Write zu and zw data (vertical axes)
3767             ALLOCATE( netcdf_data(0:nz+1) )
3768
3769             netcdf_data(0:nz+1) = zu(nzb:nzt+1)
3770             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zu_xz(av), &
3771                                     netcdf_data, start = (/ 1 /),    &
3772                                     count = (/ nz+2 /) )
3773             CALL netcdf_handle_error( 'netcdf_define_header', 166 )
3774
3775             netcdf_data(0:nz+1) = zw(nzb:nzt+1)
3776             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zw_xz(av), &
3777                                     netcdf_data, start = (/ 1 /),    &
3778                                     count = (/ nz+2 /) )
3779             CALL netcdf_handle_error( 'netcdf_define_header', 167 )
3780
3781!
3782!--          Write zs data
3783             IF ( land_surface )  THEN
3784                netcdf_data(0:nzs-1) = - zs(nzb_soil:nzt_soil)
3785                nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zs_xz(av), &
3786                                        netcdf_data(0:nzs), start = (/ 1 /),    &
3787                                        count = (/ nzt_soil-nzb_soil+1 /) )
3788               CALL netcdf_handle_error( 'netcdf_define_header', 548 )
3789             ENDIF
3790
3791             DEALLOCATE( netcdf_data )
3792!
3793!--          Write UTM coordinates
3794             IF ( init_model%rotation_angle == 0.0_wp )  THEN
3795!
3796!--             1D in case of no rotation
3797                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
3798!
3799!--             x coordinates
3800                ALLOCATE( netcdf_data(0:nx) )
3801                DO  k = 0, 2
3802!
3803!--                Scalar grid points
3804                   IF ( k == 0 )  THEN
3805                      shift_x = 0.5
3806!
3807!--                u grid points
3808                   ELSEIF ( k == 1 )  THEN
3809                      shift_x = 0.0
3810!
3811!--                v grid points
3812                   ELSEIF ( k == 2 )  THEN
3813                      shift_x = 0.5
3814                   ENDIF
3815
3816                   DO  i = 0, nx
3817                     netcdf_data(i) = init_model%origin_x            &
3818                                    + cos_ra * ( i + shift_x ) * dx
3819                   ENDDO
3820
3821                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_eutm_xz(k,av),&
3822                                           netcdf_data, start = (/ 1 /),   &
3823                                           count = (/ nx+1 /) )
3824                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
3825
3826                ENDDO
3827                DEALLOCATE( netcdf_data )
3828!
3829!--             y coordinates
3830                ALLOCATE( netcdf_data(1:ns) )
3831                DO  k = 0, 2
3832!
3833!--                Scalar grid points
3834                   IF ( k == 0 )  THEN
3835                      shift_y = 0.5
3836!
3837!--                u grid points
3838                   ELSEIF ( k == 1 )  THEN
3839                      shift_y = 0.5
3840!
3841!--                v grid points
3842                   ELSEIF ( k == 2 )  THEN
3843                      shift_y = 0.0
3844                   ENDIF
3845
3846                   DO  i = 1, ns
3847                      IF( section(i,2) == -1 )  THEN
3848                         netcdf_data(i) = -1.0_wp  ! section averaged along y
3849                      ELSE
3850                         netcdf_data(i) = init_model%origin_y &
3851                                     + cos_ra * ( section(i,2) + shift_y ) * dy
3852                      ENDIF
3853                   ENDDO
3854
3855                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_nutm_xz(k,av),&
3856                                           netcdf_data, start = (/ 1 /),   &
3857                                           count = (/ ns /) )
3858                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
3859
3860                ENDDO
3861                DEALLOCATE( netcdf_data )
3862
3863             ELSE
3864!
3865!--             2D in case of rotation
3866                ALLOCATE( netcdf_data_2d(0:nx,1:ns) )
3867                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
3868                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
3869
3870                DO  k = 0, 2
3871!
3872!--                Scalar grid points
3873                   IF ( k == 0 )  THEN
3874                      shift_x = 0.5 ; shift_y = 0.5
3875!
3876!--                u grid points
3877                   ELSEIF ( k == 1 )  THEN
3878                      shift_x = 0.0 ; shift_y = 0.5
3879!
3880!--                v grid points
3881                   ELSEIF ( k == 2 )  THEN
3882                      shift_x = 0.5 ; shift_y = 0.0
3883                   ENDIF
3884
3885                   DO  j = 1, ns
3886                      IF( section(j,2) == -1 )  THEN
3887                         netcdf_data_2d(:,j) = -1.0_wp  ! section averaged along y
3888                      ELSE
3889                         DO  i = 0, nx
3890                            netcdf_data_2d(i,j) = init_model%origin_x          &
3891                                    + cos_ra * ( i + shift_x ) * dx            &
3892                                    + sin_ra * ( section(j,2) + shift_y ) * dy
3893                         ENDDO
3894                      ENDIF
3895                   ENDDO
3896
3897                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_eutm_xz(k,av),  &
3898                                           netcdf_data_2d, start = (/ 1, 1 /),   &
3899                                           count = (/ nx+1, ns /) )
3900                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
3901
3902                   DO  j = 1, ns
3903                      IF( section(j,2) == -1 )  THEN
3904                         netcdf_data_2d(:,j) = -1.0_wp  ! section averaged along y
3905                      ELSE
3906                         DO  i = 0, nx
3907                            netcdf_data_2d(i,j) = init_model%origin_y          &
3908                                    - sin_ra * ( i + shift_x ) * dx            &
3909                                    + cos_ra * ( section(j,2) + shift_y ) * dy
3910                         ENDDO
3911                      ENDIF
3912                   ENDDO
3913
3914                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_nutm_xz(k,av),  &
3915                                           netcdf_data_2d, start = (/ 1, 1 /),   &
3916                                           count = (/ nx+1, ns /) )
3917                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
3918
3919                ENDDO
3920                DEALLOCATE( netcdf_data_2d )
3921             ENDIF
3922!
3923!--          Write lon and lat data
3924             ALLOCATE( lat(0:nx,1:ns) )
3925             ALLOCATE( lon(0:nx,1:ns) )
3926             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
3927             sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
3928
3929             DO  k = 0, 2
3930!
3931!--             Scalar grid points
3932                IF ( k == 0 )  THEN
3933                   shift_x = 0.5 ; shift_y = 0.5
3934!
3935!--             u grid points
3936                ELSEIF ( k == 1 )  THEN
3937                   shift_x = 0.0 ; shift_y = 0.5
3938!
3939!--             v grid points
3940                ELSEIF ( k == 2 )  THEN
3941                   shift_x = 0.5 ; shift_y = 0.0
3942                ENDIF
3943
3944                DO  j = 1, ns
3945                   IF( section(j,2) == -1 )  THEN
3946                      lat(:,j) = -90.0_wp  ! section averaged along y
3947                      lon(:,j) = -180.0_wp  ! section averaged along y
3948                   ELSE
3949                      DO  i = 0, nx
3950                         eutm = init_model%origin_x            &
3951                              + cos_ra * ( i + shift_x ) * dx  &
3952                              + sin_ra * ( section(j,2) + shift_y ) * dy
3953                         nutm = init_model%origin_y            &
3954                              - sin_ra * ( i + shift_x ) * dx  &
3955                              + cos_ra * ( section(j,2) + shift_y ) * dy
3956
3957                         CALL  convert_utm_to_geographic( crs_list,          &
3958                                                          eutm, nutm,        &
3959                                                          lon(i,j), lat(i,j) )
3960                      ENDDO
3961                   ENDIF
3962                ENDDO
3963
3964                nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_lon_xz(k,av), &
3965                                     lon, start = (/ 1, 1 /),       &
3966                                     count = (/ nx+1, ns /) )
3967                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
3968
3969                nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_lat_xz(k,av), &
3970                                     lat, start = (/ 1, 1 /),       &
3971                                     count = (/ nx+1, ns /) )
3972                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
3973             ENDDO
3974
3975             DEALLOCATE( lat )
3976             DEALLOCATE( lon )
3977
3978          ENDIF
3979
3980
3981       CASE ( 'xz_ext' )
3982
3983!
3984!--       Get the list of variables and compare with the actual run.
3985!--       First var_list_old has to be reset, since GET_ATT does not assign
3986!--       trailing blanks.
3987          var_list_old = ' '
3988          nc_stat = NF90_GET_ATT( id_set_xz(av), NF90_GLOBAL, 'VAR_LIST', &
3989                                  var_list_old )
3990          CALL netcdf_handle_error( 'netcdf_define_header', 168 )
3991
3992          var_list = ';'
3993          i = 1
3994          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
3995             IF ( INDEX( do2d(av,i), 'xz' ) /= 0 )  THEN
3996                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
3997             ENDIF
3998             i = i + 1
3999          ENDDO
4000
4001          IF ( av == 0 )  THEN
4002             var = '(xz)'
4003          ELSE
4004             var = '(xz_av)'
4005          ENDIF
4006
4007          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
4008             message_string = 'netCDF file for cross-sections ' //           &
4009                              TRIM( var ) // ' from previous run found,' //  &
4010                              '&but this file cannot be extended due to' //  &
4011                              ' variable mismatch.' //                       &
4012                              '&New file is created instead.'
4013             CALL message( 'define_netcdf_header', 'PA0249', 0, 1, 0, 6, 0 )
4014             extend = .FALSE.
4015             RETURN
4016          ENDIF
4017
4018!
4019!--       Calculate the number of current sections
4020          ns = 1
4021          DO WHILE ( section(ns,2) /= -9999  .AND.  ns <= 100 )
4022             ns = ns + 1
4023          ENDDO
4024          ns = ns - 1
4025
4026!
4027!--       Get and compare the number of vertical cross sections
4028          nc_stat = NF90_INQ_VARID( id_set_xz(av), 'y_xz', id_var_y_xz(av) )
4029          CALL netcdf_handle_error( 'netcdf_define_header', 169 )
4030
4031          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_y_xz(av), &
4032                                           dimids = id_dim_y_xz_old )
4033          CALL netcdf_handle_error( 'netcdf_define_header', 170 )
4034          id_dim_y_xz(av) = id_dim_y_xz_old(1)
4035
4036          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_y_xz(av), &
4037                                            len = ns_old )
4038          CALL netcdf_handle_error( 'netcdf_define_header', 171 )
4039
4040          IF ( ns /= ns_old )  THEN
4041             message_string = 'netCDF file for cross-sections ' //          &
4042                              TRIM( var ) // ' from previous run found,' // &
4043                              '&but this file cannot be extended due to' // &
4044                              ' mismatch in number of' //                   &
4045                              ' cross sections.' //                         &
4046                              '&New file is created instead.'
4047             CALL message( 'define_netcdf_header', 'PA0250', 0, 1, 0, 6, 0 )
4048             extend = .FALSE.
4049             RETURN
4050          ENDIF
4051
4052!
4053!--       Get and compare the heights of the cross sections
4054          ALLOCATE( netcdf_data(1:ns_old) )
4055
4056          nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_y_xz(av), netcdf_data )
4057          CALL netcdf_handle_error( 'netcdf_define_header', 172 )
4058
4059          DO  i = 1, ns
4060             IF ( section(i,2) /= -1 )  THEN
4061                IF ( ( ( section(i,2) + 0.5 ) * dy ) /= netcdf_data(i) )  THEN
4062                   message_string = 'netCDF file for cross-sections ' //       &
4063                               TRIM( var ) // ' from previous run found,' //   &
4064                               ' but this file cannot be extended' //          &
4065                               ' due to mismatch in cross' //                  &
4066                               ' section levels.' //                           &
4067                               ' New file is created instead.'
4068                   CALL message( 'define_netcdf_header', 'PA0251',             &
4069                                                                 0, 1, 0, 6, 0 )
4070                   extend = .FALSE.
4071                   RETURN
4072                ENDIF
4073             ELSE
4074                IF ( -1.0_wp /= netcdf_data(i) )  THEN
4075                   message_string = 'netCDF file for cross-sections ' //       &
4076                               TRIM( var ) // ' from previous run found,' //   &
4077                               ' but this file cannot be extended' //          &
4078                               ' due to mismatch in cross' //                  &
4079                               ' section levels.' //                           &
4080                               ' New file is created instead.'
4081                   CALL message( 'define_netcdf_header', 'PA0251',             &
4082                                                                 0, 1, 0, 6, 0 )
4083                   extend = .FALSE.
4084                   RETURN
4085                ENDIF
4086             ENDIF
4087          ENDDO
4088
4089          DEALLOCATE( netcdf_data )
4090
4091!
4092!--       Get the id of the time coordinate (unlimited coordinate) and its
4093!--       last index on the file. The next time level is do2d..count+1.
4094!--       The current time must be larger than the last output time
4095!--       on the file.
4096          nc_stat = NF90_INQ_VARID( id_set_xz(av), 'time', id_var_time_xz(av) )
4097          CALL netcdf_handle_error( 'netcdf_define_header', 173 )
4098
4099          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_time_xz(av), &
4100                                           dimids = id_dim_time_old )
4101          CALL netcdf_handle_error( 'netcdf_define_header', 174 )
4102          id_dim_time_xz(av) = id_dim_time_old(1)
4103
4104          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_time_xz(av), &
4105                                            len = ntime_count )
4106          CALL netcdf_handle_error( 'netcdf_define_header', 175 )
4107
4108!
4109!--       For non-parallel output use the last output time level of the netcdf
4110!--       file because the time dimension is unlimited. In case of parallel
4111!--       output the variable ntime_count could get the value of 9*10E36 because
4112!--       the time dimension is limited.
4113          IF ( netcdf_data_format < 5 ) do2d_xz_time_count(av) = ntime_count
4114
4115          nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_time_xz(av),           &
4116                                  last_time_coordinate,                        &
4117                                  start = (/ do2d_xz_time_count(av) /),        &
4118                                  count = (/ 1 /) )
4119          CALL netcdf_handle_error( 'netcdf_define_header', 176 )
4120
4121          IF ( last_time_coordinate(1) >= simulated_time )  THEN
4122             message_string = 'netCDF file for cross sections ' //             &
4123                              TRIM( var ) // ' from previous run found,' //    &
4124                              '&but this file cannot be extended becaus' //    &
4125                              'e the current output time' //                   &
4126                              '&is less or equal than the last output t' //    &
4127                              'ime on this file.' //                           &
4128                              '&New file is created instead.'
4129             CALL message( 'define_netcdf_header', 'PA0252', 0, 1, 0, 6, 0 )
4130             do2d_xz_time_count(av) = 0
4131             extend = .FALSE.
4132             RETURN
4133          ENDIF
4134
4135          IF ( netcdf_data_format > 4 )  THEN
4136!
4137!--          Check if the needed number of output time levels is increased
4138!--          compared to the number of time levels in the existing file.
4139             IF ( ntdim_2d_xz(av) > ntime_count )  THEN
4140                message_string = 'netCDF file for cross sections ' // &
4141                                 TRIM( var ) // ' from previous run found,' // &
4142                                 '&but this file cannot be extended becaus' // &
4143                                 'e the number of output time levels has b' // &
4144                                 'een increased compared to the previous s' // &
4145                                 'imulation.' //                               &
4146                                 '&New file is created instead.'
4147                CALL message( 'define_netcdf_header', 'PA0390', 0, 1, 0, 6, 0 )
4148                do2d_xz_time_count(av) = 0
4149                extend = .FALSE.
4150!
4151!--             Recalculate the needed time levels for the new file.
4152                IF ( av == 0 )  THEN
4153                   ntdim_2d_xz(0) = CEILING(                            &
4154                           ( end_time - MAX( skip_time_do2d_xz,         &
4155                                             simulated_time_at_begin )  &
4156                           ) / dt_do2d_xz )
4157                   IF ( do2d_at_begin )  ntdim_2d_xz(0) = ntdim_2d_xz(0) + 1
4158                ELSE
4159                   ntdim_2d_xz(1) = CEILING(                            &
4160                           ( end_time - MAX( skip_time_data_output_av,  &
4161                                             simulated_time_at_begin )  &
4162                           ) / dt_data_output_av )
4163                ENDIF
4164                RETURN
4165             ENDIF
4166          ENDIF
4167
4168!
4169!--       Dataset seems to be extendable.
4170!--       Now get the variable ids.
4171          i = 1
4172          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
4173             IF ( INDEX( do2d(av,i), 'xz' ) /= 0 )  THEN
4174                nc_stat = NF90_INQ_VARID( id_set_xz(av), do2d(av,i), &
4175                                          id_var_do2d(av,i) )
4176                CALL netcdf_handle_error( 'netcdf_define_header', 177 )
4177#if defined( __netcdf4_parallel )
4178!
4179!--             Set independent io operations for parallel io. Collective io
4180!--             is only allowed in case of a 1d-decomposition along x, because
4181!--             otherwise, not all PEs have output data.
4182                IF ( netcdf_data_format > 4 )  THEN
4183                   IF ( npey == 1 )  THEN
4184                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
4185                                                     id_var_do2d(av,i), &
4186                                                     NF90_COLLECTIVE )
4187                   ELSE
4188!
4189!--                   Test simulations showed that the output of cross sections
4190!--                   by all PEs in data_output_2d using NF90_COLLECTIVE is
4191!--                   faster than the output by the first row of PEs in
4192!--                   x-direction using NF90_INDEPENDENT.
4193                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
4194                                                     id_var_do2d(av,i), &
4195                                                     NF90_COLLECTIVE )
4196!                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
4197!                                                     id_var_do2d(av,i), &
4198!                                                     NF90_INDEPENDENT )
4199                   ENDIF
4200                   CALL netcdf_handle_error( 'netcdf_define_header', 455 )
4201                ENDIF
4202#endif
4203             ENDIF
4204             i = i + 1
4205          ENDDO
4206
4207!
4208!--       Update the title attribute on file
4209!--       In order to avoid 'data mode' errors if updated attributes are larger
4210!--       than their original size, NF90_PUT_ATT is called in 'define mode'
4211!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
4212!--       performance loss due to data copying; an alternative strategy would be
4213!--       to ensure equal attribute size in a job chain. Maybe revise later.
4214          IF ( av == 0 )  THEN
4215             time_average_text = ' '
4216          ELSE
4217             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
4218                                                            averaging_interval
4219          ENDIF
4220          nc_stat = NF90_REDEF( id_set_xz(av) )
4221          CALL netcdf_handle_error( 'netcdf_define_header', 433 )
4222          nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'title',         &
4223                                  TRIM( run_description_header ) //            &
4224                                  TRIM( time_average_text ) )
4225          CALL netcdf_handle_error( 'netcdf_define_header', 178 )
4226          nc_stat = NF90_ENDDEF( id_set_xz(av) )
4227          CALL netcdf_handle_error( 'netcdf_define_header', 434 )
4228          message_string = 'netCDF file for cross-sections ' //                &
4229                            TRIM( var ) // ' from previous run found.' //      &
4230                           '&This file will be extended.'
4231          CALL message( 'define_netcdf_header', 'PA0253', 0, 0, 0, 6, 0 )
4232
4233
4234       CASE ( 'yz_new' )
4235
4236!
4237!--       Define some global attributes of the dataset
4238          IF ( av == 0 )  THEN
4239             CALL netcdf_create_global_atts( id_set_yz(av), 'yz', TRIM( run_description_header ), 179 )
4240             time_average_text = ' '
4241          ELSE
4242             CALL netcdf_create_global_atts( id_set_yz(av), 'yz_av', TRIM( run_description_header ), 179 )
4243             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
4244             nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'time_avg',   &
4245                                     TRIM( time_average_text ) )
4246             CALL netcdf_handle_error( 'netcdf_define_header', 180 )
4247          ENDIF
4248
4249!
4250!--       Define time coordinate for yz sections.
4251!--       For parallel output the time dimensions has to be limited, otherwise
4252!--       the performance drops significantly.
4253          IF ( netcdf_data_format < 5 )  THEN
4254             CALL netcdf_create_dim( id_set_yz(av), 'time', NF90_UNLIMITED,    &
4255                                     id_dim_time_yz(av), 181 )
4256          ELSE
4257             CALL netcdf_create_dim( id_set_yz(av), 'time', ntdim_2d_yz(av),   &
4258                                     id_dim_time_yz(av), 526 )
4259          ENDIF
4260
4261          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_time_yz(av) /),     &
4262                                  'time', NF90_DOUBLE, id_var_time_yz(av),     &
4263                                  'seconds', 'time', 182, 183, 000 )
4264          CALL netcdf_create_att( id_set_yz(av), id_var_time_yz(av), 'standard_name', 'time', 000)
4265          CALL netcdf_create_att( id_set_yz(av), id_var_time_yz(av), 'axis', 'T', 000)
4266!
4267!--       Define the spatial dimensions and coordinates for yz-sections.
4268!--       First, determine the number of vertical sections to be written.
4269          IF ( section(1,3) == -9999 )  THEN
4270             RETURN
4271          ELSE
4272             ns = 1
4273             DO WHILE ( section(ns,3) /= -9999  .AND.  ns <= 100 )
4274                ns = ns + 1
4275             ENDDO
4276             ns = ns - 1
4277          ENDIF
4278
4279!
4280!--       Define x axis (for scalar position)
4281          CALL netcdf_create_dim( id_set_yz(av), 'x_yz', ns, id_dim_x_yz(av),  &
4282                                  184 )
4283          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_x_yz(av) /),        &
4284                                  'x_yz', NF90_DOUBLE, id_var_x_yz(av),        &
4285                                  'meters', '', 185, 186, 000 )
4286!
4287!--       Define x axis (for u position)
4288          CALL netcdf_create_dim( id_set_yz(av), 'xu_yz', ns,                  &
4289                                  id_dim_xu_yz(av), 377 )
4290          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_xu_yz(av) /),       &
4291                                  'xu_yz', NF90_DOUBLE, id_var_xu_yz(av),      &
4292                                  'meters', '', 378, 379, 000 )
4293!
4294!--       Define a variable to store the layer indices of the vertical cross
4295!--       sections
4296          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_x_yz(av) /),        &
4297                                  'ind_x_yz', NF90_DOUBLE,                     &
4298                                  id_var_ind_x_yz(av), 'gridpoints', '', 187,  &
4299                                  188, 000 )
4300!
4301!--       Define y-axis (for scalar position)
4302          CALL netcdf_create_dim( id_set_yz(av), 'y', ny+1, id_dim_y_yz(av),   &
4303                                  189 )
4304          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_y_yz(av) /), 'y',   &
4305                                  NF90_DOUBLE, id_var_y_yz(av), 'meters', '',  &
4306                                  190, 191, 000 )
4307!
4308!--       Define y-axis (for v position)
4309          CALL netcdf_create_dim( id_set_yz(av), 'yv', ny+1, id_dim_yv_yz(av), &
4310                                  380 )
4311          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_yv_yz(av) /), 'yv', &
4312                                  NF90_DOUBLE, id_var_yv_yz(av), 'meters', '', &
4313                                  381, 382, 000 )
4314!
4315!--       Define the two z-axes (zu and zw)
4316          CALL netcdf_create_dim( id_set_yz(av), 'zu', nz+2, id_dim_zu_yz(av), &
4317                                  192 )
4318          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_zu_yz(av) /), 'zu', &
4319                                  NF90_DOUBLE, id_var_zu_yz(av), 'meters', '', &
4320                                  193, 194, 000 )
4321
4322          CALL netcdf_create_dim( id_set_yz(av), 'zw', nz+2, id_dim_zw_yz(av), &
4323                                  195 )
4324          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_zw_yz(av) /), 'zw', &
4325                                  NF90_DOUBLE, id_var_zw_yz(av), 'meters', '', &
4326                                  196, 197, 000 )
4327!
4328!--       Define UTM and geographic coordinates
4329          CALL define_geo_coordinates( id_set_yz(av),         &
4330                  (/ id_dim_x_yz(av), id_dim_xu_yz(av) /),    &
4331                  (/ id_dim_y_yz(av), id_dim_yv_yz(av) /),    &
4332                  id_var_eutm_yz(:,av), id_var_nutm_yz(:,av), &
4333                  id_var_lat_yz(:,av), id_var_lon_yz(:,av)    )
4334!
4335!--       Define coordinate-reference system
4336          CALL netcdf_create_crs( id_set_yz(av), 000 )
4337
4338          IF ( land_surface )  THEN
4339
4340             CALL netcdf_create_dim( id_set_yz(av), 'zs', nzs,                 &
4341                                     id_dim_zs_yz(av), 545 )
4342             CALL netcdf_create_var( id_set_yz(av), (/ id_dim_zs_yz(av) /),    &
4343                                     'zs', NF90_DOUBLE, id_var_zs_yz(av),      &
4344                                     'meters', '', 546, 547, 000 )
4345
4346          ENDIF
4347
4348!
4349!--       Define the variables
4350          var_list = ';'
4351          i = 1
4352
4353          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
4354
4355             IF ( INDEX( do2d(av,i), 'yz' ) /= 0 )  THEN
4356
4357!
4358!--             Check for the grid
4359                found = .FALSE.
4360                SELECT CASE ( do2d(av,i) )
4361!
4362!--                Most variables are defined on the zu grid
4363                   CASE ( 'e_yz', 'nc_yz', 'nr_yz', 'p_yz', 'pc_yz',           &
4364                          'pr_yz','prr_yz', 'q_yz', 'qc_yz', 'ql_yz',          &
4365                          'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz',  &
4366                          's_yz',                                              &
4367                          'theta_yz', 'thetal_yz', 'thetav_yz', 'ti_yz' )
4368
4369                      grid_x = 'x'
4370                      grid_y = 'y'
4371                      grid_z = 'zu'
4372!
4373!--                u grid
4374                   CASE ( 'u_yz' )
4375
4376                      grid_x = 'xu'
4377                      grid_y = 'y'
4378                      grid_z = 'zu'
4379!
4380!--                v grid
4381                   CASE ( 'v_yz' )
4382
4383                      grid_x = 'x'
4384                      grid_y = 'yv'
4385                      grid_z = 'zu'
4386!
4387!--                w grid
4388                   CASE ( 'w_yz' )
4389
4390                      grid_x = 'x'
4391                      grid_y = 'y'
4392                      grid_z = 'zw'
4393
4394
4395                   CASE DEFAULT
4396!
4397!--                   Check for land surface quantities
4398                      IF ( land_surface )  THEN
4399                         CALL lsm_define_netcdf_grid( do2d(av,i), found,       &
4400                                                      grid_x, grid_y, grid_z )
4401                      ENDIF
4402
4403                      IF ( .NOT. found )  THEN
4404                         CALL tcm_define_netcdf_grid( do2d(av,i), found,       &
4405                                                      grid_x, grid_y, grid_z )
4406                      ENDIF
4407
4408!
4409!--                   Check for ocean quantities
4410                      IF ( .NOT. found  .AND.  ocean_mode )  THEN
4411                         CALL ocean_define_netcdf_grid( do2d(av,i), found,     &
4412                                                       grid_x, grid_y, grid_z )
4413                      ENDIF
4414!
4415!--                   Check for radiation quantities
4416                      IF ( .NOT. found  .AND.  radiation )  THEN
4417                         CALL radiation_define_netcdf_grid( do2d(av,i), found, &
4418                                                            grid_x, grid_y,    &
4419                                                            grid_z )
4420                      ENDIF
4421!
4422!--                   Check for SALSA quantities
4423                      IF ( .NOT. found  .AND.  salsa )  THEN
4424                         CALL salsa_define_netcdf_grid( do2d(av,i), found,     &
4425                                                        grid_x, grid_y, grid_z )
4426                      ENDIF
4427!
4428!--                   Check for gust module quantities
4429                      IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
4430                         CALL gust_define_netcdf_grid( do2d(av,i), found,      &
4431                                                       grid_x, grid_y, grid_z )
4432                      ENDIF
4433
4434!
4435!--                   Check for chemistry quantities
4436                      IF ( .NOT. found  .AND.  air_chemistry )  THEN
4437                         CALL chem_define_netcdf_grid( do2d(av,i), found,      &
4438                                                       grid_x, grid_y,         &
4439                                                       grid_z )
4440                      ENDIF
4441
4442                      IF ( .NOT. found )                                       &
4443                         CALL doq_define_netcdf_grid(                          &
4444                                                    do2d(av,i), found, grid_x, &
4445                                                    grid_y, grid_z           )
4446!
4447!--                   Check for user-defined quantities
4448                      IF ( .NOT. found  .AND.  user_module_enabled )  THEN
4449                         CALL user_define_netcdf_grid( do2d(av,i), found,      &
4450                                                       grid_x, grid_y, grid_z )
4451                      ENDIF
4452
4453                      IF ( .NOT. found )  THEN
4454                         WRITE ( message_string, * ) 'no grid defined for',    &
4455                                                ' variable ', TRIM( do2d(av,i) )
4456                         CALL message( 'define_netcdf_header', 'PA0244',       &
4457                                       0, 1, 0, 6, 0 )
4458                      ENDIF
4459
4460                END SELECT
4461
4462!
4463!--             Select the respective dimension ids
4464                IF ( grid_x == 'x' )  THEN
4465                   id_x = id_dim_x_yz(av)
4466                ELSEIF ( grid_x == 'xu' )  THEN
4467                   id_x = id_dim_xu_yz(av)
4468                ENDIF
4469
4470                IF ( grid_y == 'y' )  THEN
4471                   id_y = id_dim_y_yz(av)
4472                ELSEIF ( grid_y == 'yv' )  THEN
4473                   id_y = id_dim_yv_yz(av)
4474                ENDIF
4475
4476                IF ( grid_z == 'zu' )  THEN
4477                   id_z = id_dim_zu_yz(av)
4478                ELSEIF ( grid_z == 'zw' )  THEN
4479                   id_z = id_dim_zw_yz(av)
4480                ELSEIF ( grid_z == 'zs' )  THEN
4481                   id_z = id_dim_zs_yz(av)
4482                ENDIF
4483
4484!
4485!--             Define the grid
4486                CALL netcdf_create_var( id_set_yz(av),  (/ id_x, id_y, id_z,   &
4487                                        id_dim_time_yz(av) /), do2d(av,i),     &
4488                                        nc_precision(3), id_var_do2d(av,i),    &
4489                                        TRIM( do2d_unit(av,i) ), do2d(av,i),   &
4490                                        198, 199, 356, .TRUE. )
4491
4492#if defined( __netcdf4_parallel )
4493                IF ( netcdf_data_format > 4 )  THEN
4494!
4495!--                Set no fill for every variable to increase performance.
4496                   nc_stat = NF90_DEF_VAR_FILL( id_set_yz(av),     &
4497                                                id_var_do2d(av,i), &
4498                                                NF90_NOFILL, 0 )
4499                   CALL netcdf_handle_error( 'netcdf_define_header', 535 )
4500!
4501!--                Set independent io operations for parallel io. Collective io
4502!--                is only allowed in case of a 1d-decomposition along y,
4503!--                because otherwise, not all PEs have output data.
4504                   IF ( npex == 1 )  THEN
4505                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av),     &
4506                                                     id_var_do2d(av,i), &
4507                                                     NF90_COLLECTIVE )
4508                   ELSE
4509!
4510!--                   Test simulations showed that the output of cross sections
4511!--                   by all PEs in data_output_2d using NF90_COLLECTIVE is
4512!--                   faster than the output by the first row of PEs in
4513!--                   y-direction using NF90_INDEPENDENT.
4514                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av),     &
4515                                                     id_var_do2d(av,i), &
4516                                                     NF90_COLLECTIVE )
4517!                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av),     &
4518!                                                     id_var_do2d(av,i), &
4519!                                                     NF90_INDEPENDENT )
4520                   ENDIF
4521                   CALL netcdf_handle_error( 'netcdf_define_header', 450 )
4522                ENDIF
4523#endif
4524                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
4525
4526             ENDIF
4527
4528             i = i + 1
4529
4530          ENDDO
4531
4532!
4533!--       No arrays to output. Close the netcdf file and return.
4534          IF ( i == 1 )  RETURN
4535
4536!
4537!--       Write the list of variables as global attribute (this is used by
4538!--       restart runs and by combine_plot_fields)
4539          nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'VAR_LIST', &
4540                                  var_list )
4541          CALL netcdf_handle_error( 'netcdf_define_header', 200 )
4542
4543!
4544!--       Set general no fill, otherwise the performance drops significantly for
4545!--       parallel output.
4546          nc_stat = NF90_SET_FILL( id_set_yz(av), NF90_NOFILL, oldmode )
4547          CALL netcdf_handle_error( 'netcdf_define_header', 531 )
4548
4549!
4550!--       Leave netCDF define mode
4551          nc_stat = NF90_ENDDEF( id_set_yz(av) )
4552          CALL netcdf_handle_error( 'netcdf_define_header', 201 )
4553
4554!
4555!--       These data are only written by PE0 for parallel output to increase
4556!--       the performance.
4557          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
4558
4559!
4560!--          Write axis data: x_yz, y, zu, zw
4561             ALLOCATE( netcdf_data(1:ns) )
4562
4563!
4564!--          Write x_yz data (shifted by +dx/2)
4565             DO  i = 1, ns
4566                IF( section(i,3) == -1 )  THEN
4567                   netcdf_data(i) = -1.0_wp  ! section averaged along x
4568                ELSE
4569                   netcdf_data(i) = ( section(i,3) + 0.5_wp ) * dx
4570                ENDIF
4571             ENDDO
4572             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_x_yz(av), &
4573                                     netcdf_data, start = (/ 1 /),   &
4574                                     count = (/ ns /) )
4575             CALL netcdf_handle_error( 'netcdf_define_header', 202 )
4576
4577!
4578!--          Write x_yz data (xu grid)
4579             DO  i = 1, ns
4580                IF( section(i,3) == -1 )  THEN
4581                   netcdf_data(i) = -1.0_wp  ! section averaged along x
4582                ELSE
4583                   netcdf_data(i) = section(i,3) * dx
4584                ENDIF
4585             ENDDO
4586             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_xu_yz(av), &
4587                                     netcdf_data, start = (/ 1 /),    &
4588                                     count = (/ ns /) )
4589             CALL netcdf_handle_error( 'netcdf_define_header', 383 )
4590
4591!
4592!--          Write gridpoint number data
4593             netcdf_data(1:ns) = section(1:ns,3)
4594             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_ind_x_yz(av), &
4595                                     netcdf_data, start = (/ 1 /),       &
4596                                     count = (/ ns /) )
4597             CALL netcdf_handle_error( 'netcdf_define_header', 203 )
4598
4599             DEALLOCATE( netcdf_data )
4600
4601!
4602!--          Write data for y (shifted by +dy/2) and yv axis
4603             ALLOCATE( netcdf_data(0:ny) )
4604
4605             DO  j = 0, ny
4606                netcdf_data(j) = ( j + 0.5_wp ) * dy
4607             ENDDO
4608
4609             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_y_yz(av), &
4610                                     netcdf_data, start = (/ 1 /),   &
4611                                     count = (/ ny+1 /) )
4612             CALL netcdf_handle_error( 'netcdf_define_header', 204 )
4613
4614             DO  j = 0, ny
4615                netcdf_data(j) = j * dy
4616             ENDDO
4617
4618             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_yv_yz(av), &
4619                                     netcdf_data, start = (/ 1 /),    &
4620                                     count = (/ ny+1 /) )
4621             CALL netcdf_handle_error( 'netcdf_define_header', 384 )
4622
4623             DEALLOCATE( netcdf_data )
4624
4625!
4626!--          Write zu and zw data (vertical axes)
4627             ALLOCATE( netcdf_data(0:nz+1) )
4628
4629             netcdf_data(0:nz+1) = zu(nzb:nzt+1)
4630             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zu_yz(av), &
4631                                     netcdf_data, start = (/ 1 /),    &
4632                                     count = (/ nz+2 /) )
4633             CALL netcdf_handle_error( 'netcdf_define_header', 205 )
4634
4635             netcdf_data(0:nz+1) = zw(nzb:nzt+1)
4636             nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zw_yz(av), &
4637                                     netcdf_data, start = (/ 1 /),    &
4638                                     count = (/ nz+2 /) )
4639             CALL netcdf_handle_error( 'netcdf_define_header', 206 )
4640
4641             DEALLOCATE( netcdf_data )
4642!
4643!--          Write UTM coordinates
4644             IF ( init_model%rotation_angle == 0.0_wp )  THEN
4645!
4646!--             1D in case of no rotation
4647                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
4648!
4649!--             x coordinates
4650                ALLOCATE( netcdf_data(1:ns) )
4651                DO  k = 0, 2
4652!
4653!--                Scalar grid points
4654                   IF ( k == 0 )  THEN
4655                      shift_x = 0.5
4656!
4657!--                u grid points
4658                   ELSEIF ( k == 1 )  THEN
4659                      shift_x = 0.0
4660!
4661!--                v grid points
4662                   ELSEIF ( k == 2 )  THEN
4663                      shift_x = 0.5
4664                   ENDIF
4665
4666                   DO  i = 1, ns
4667                      IF( section(i,3) == -1 )  THEN
4668                         netcdf_data(i) = -1.0_wp  ! section averaged along x
4669                      ELSE
4670                         netcdf_data(i) = init_model%origin_x &
4671                                     + cos_ra * ( section(i,3) + shift_x ) * dx
4672                      ENDIF
4673                   ENDDO
4674
4675                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_eutm_yz(k,av),&
4676                                           netcdf_data, start = (/ 1 /),   &
4677                                           count = (/ ns /) )
4678                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
4679
4680                ENDDO
4681                DEALLOCATE( netcdf_data )
4682!
4683!--             y coordinates
4684                ALLOCATE( netcdf_data(0:ny) )
4685                DO  k = 0, 2
4686!
4687!--                Scalar grid points
4688                   IF ( k == 0 )  THEN
4689                      shift_y = 0.5
4690!
4691!--                u grid points
4692                   ELSEIF ( k == 1 )  THEN
4693                      shift_y = 0.5
4694!
4695!--                v grid points
4696                   ELSEIF ( k == 2 )  THEN
4697                      shift_y = 0.0
4698                   ENDIF
4699
4700                   DO  i = 0, ny
4701                     netcdf_data(i) = init_model%origin_y            &
4702                                    + cos_ra * ( i + shift_y ) * dy
4703                   ENDDO
4704
4705                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_nutm_yz(k,av),&
4706                                           netcdf_data, start = (/ 1 /),   &
4707                                           count = (/ ny+1 /) )
4708                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
4709
4710                ENDDO
4711                DEALLOCATE( netcdf_data )
4712
4713             ELSE
4714!
4715!--             2D in case of rotation
4716                ALLOCATE( netcdf_data_2d(1:ns,0:ny) )
4717                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
4718                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
4719
4720                DO  k = 0, 2
4721!
4722!--                Scalar grid points
4723                   IF ( k == 0 )  THEN
4724                      shift_x = 0.5 ; shift_y = 0.5
4725!
4726!--                u grid points
4727                   ELSEIF ( k == 1 )  THEN
4728                      shift_x = 0.0 ; shift_y = 0.5
4729!
4730!--                v grid points
4731                   ELSEIF ( k == 2 )  THEN
4732                      shift_x = 0.5 ; shift_y = 0.0
4733                   ENDIF
4734
4735                   DO  j = 0, ny
4736                      DO  i = 1, ns
4737                         IF( section(i,3) == -1 )  THEN
4738                            netcdf_data_2d(i,:) = -1.0_wp !section averaged along x
4739                         ELSE
4740                            netcdf_data_2d(i,j) = init_model%origin_x          &
4741                                    + cos_ra * ( section(i,3) + shift_x ) * dx &
4742                                    + sin_ra * ( j + shift_y ) * dy
4743                         ENDIF
4744                      ENDDO
4745                   ENDDO
4746
4747                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_eutm_yz(k,av),  &
4748                                           netcdf_data_2d, start = (/ 1, 1 /),   &
4749                                           count = (/ ns, ny+1 /) )
4750                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
4751
4752                   DO  j = 0, ny
4753                      DO  i = 1, ns
4754                         IF( section(i,3) == -1 )  THEN
4755                            netcdf_data_2d(i,:) = -1.0_wp !section averaged along x
4756                         ELSE
4757                            netcdf_data_2d(i,j) = init_model%origin_y          &
4758                                    - sin_ra * ( section(i,3) + shift_x ) * dx &
4759                                    + cos_ra * ( j + shift_y ) * dy
4760                         ENDIF
4761                      ENDDO
4762                   ENDDO
4763
4764                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_nutm_yz(k,av),  &
4765                                           netcdf_data_2d, start = (/ 1, 1 /),   &
4766                                           count = (/ ns, ny+1 /) )
4767                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
4768
4769                ENDDO
4770                DEALLOCATE( netcdf_data_2d )
4771             ENDIF
4772!
4773!--          Write lon and lat data
4774             ALLOCATE( lat(1:ns,0:ny) )
4775             ALLOCATE( lon(1:ns,0:ny) )
4776             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
4777             sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
4778
4779             DO  k = 0, 2
4780!
4781!--             Scalar grid points
4782                IF ( k == 0 )  THEN
4783                   shift_x = 0.5 ; shift_y = 0.5
4784!
4785!--             u grid points
4786                ELSEIF ( k == 1 )  THEN
4787                   shift_x = 0.0 ; shift_y = 0.5
4788!
4789!--             v grid points
4790                ELSEIF ( k == 2 )  THEN
4791                   shift_x = 0.5 ; shift_y = 0.0
4792                ENDIF
4793
4794                DO  j = 0, ny
4795                   DO  i = 1, ns
4796                      IF( section(i,3) == -1 )  THEN
4797                         lat(i,:) = -90.0_wp   ! section averaged along x
4798                         lon(i,:) = -180.0_wp  ! section averaged along x
4799                      ELSE
4800                         eutm = init_model%origin_x            &
4801                              + cos_ra * ( section(i,3) + shift_x ) * dx  &
4802                              + sin_ra * ( j + shift_y ) * dy
4803                         nutm = init_model%origin_y            &
4804                              - sin_ra * ( section(i,3) + shift_x ) * dx  &
4805                              + cos_ra * ( j + shift_y ) * dy
4806
4807                         CALL  convert_utm_to_geographic( crs_list,          &
4808                                                          eutm, nutm,        &
4809                                                          lon(i,j), lat(i,j) )
4810                      ENDIF
4811                   ENDDO
4812                ENDDO
4813
4814                nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_lon_yz(k,av), &
4815                                     lon, start = (/ 1, 1 /),       &
4816                                     count = (/ ns, ny+1 /) )
4817                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
4818
4819                nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_lat_yz(k,av), &
4820                                     lat, start = (/ 1, 1 /),       &
4821                                     count = (/ ns, ny+1 /) )
4822                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
4823             ENDDO
4824
4825             DEALLOCATE( lat )
4826             DEALLOCATE( lon )
4827
4828          ENDIF
4829
4830
4831       CASE ( 'yz_ext' )
4832
4833!
4834!--       Get the list of variables and compare with the actual run.
4835!--       First var_list_old has to be reset, since GET_ATT does not assign
4836!--       trailing blanks.
4837          var_list_old = ' '
4838          nc_stat = NF90_GET_ATT( id_set_yz(av), NF90_GLOBAL, 'VAR_LIST', &
4839                                  var_list_old )
4840          CALL netcdf_handle_error( 'netcdf_define_header', 207 )
4841
4842          var_list = ';'
4843          i = 1
4844          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
4845             IF ( INDEX( do2d(av,i), 'yz' ) /= 0 )  THEN
4846                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
4847             ENDIF
4848             i = i + 1
4849          ENDDO
4850
4851          IF ( av == 0 )  THEN
4852             var = '(yz)'
4853          ELSE
4854             var = '(yz_av)'
4855          ENDIF
4856
4857          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
4858             message_string = 'netCDF file for cross-sections ' //           &
4859                              TRIM( var ) // ' from previous run found,' //  &
4860                              '&but this file cannot be extended due to' //  &
4861                              ' variable mismatch.' //                       &
4862                              '&New file is created instead.'
4863             CALL message( 'define_netcdf_header', 'PA0249', 0, 1, 0, 6, 0 )
4864             extend = .FALSE.
4865             RETURN
4866          ENDIF
4867
4868!
4869!--       Calculate the number of current sections
4870          ns = 1
4871          DO WHILE ( section(ns,3) /= -9999  .AND.  ns <= 100 )
4872             ns = ns + 1
4873          ENDDO
4874          ns = ns - 1
4875
4876!
4877!--       Get and compare the number of vertical cross sections
4878          nc_stat = NF90_INQ_VARID( id_set_yz(av), 'x_yz', id_var_x_yz(av) )
4879          CALL netcdf_handle_error( 'netcdf_define_header', 208 )
4880
4881          nc_stat = NF90_INQUIRE_VARIABLE( id_set_yz(av), id_var_x_yz(av), &
4882                                           dimids = id_dim_x_yz_old )
4883          CALL netcdf_handle_error( 'netcdf_define_header', 209 )
4884          id_dim_x_yz(av) = id_dim_x_yz_old(1)
4885
4886          nc_stat = NF90_INQUIRE_DIMENSION( id_set_yz(av), id_dim_x_yz(av), &
4887                                            len = ns_old )
4888          CALL netcdf_handle_error( 'netcdf_define_header', 210 )
4889
4890          IF ( ns /= ns_old )  THEN
4891             message_string = 'netCDF file for cross-sections ' //          &
4892                              TRIM( var ) // ' from previous run found,' // &
4893                              '&but this file cannot be extended due to' // &
4894                              ' mismatch in number of' //                   &
4895                              ' cross sections.' //                         &
4896                              '&New file is created instead.'
4897             CALL message( 'define_netcdf_header', 'PA0250', 0, 1, 0, 6, 0 )
4898             extend = .FALSE.
4899             RETURN
4900          ENDIF
4901
4902!
4903!--       Get and compare the heights of the cross sections
4904          ALLOCATE( netcdf_data(1:ns_old) )
4905
4906          nc_stat = NF90_GET_VAR( id_set_yz(av), id_var_x_yz(av), netcdf_data )
4907          CALL netcdf_handle_error( 'netcdf_define_header', 211 )
4908
4909          DO  i = 1, ns
4910             IF ( section(i,3) /= -1 )  THEN
4911                IF ( ( ( section(i,3) + 0.5 ) * dx ) /= netcdf_data(i) )  THEN
4912                   message_string = 'netCDF file for cross-sections ' //       &
4913                              TRIM( var ) // ' from previous run found,' //    &
4914                              ' but this file cannot be extended' //           &
4915                              ' due to mismatch in cross' //                   &
4916                              ' section levels.' //                            &
4917                              ' New file is created instead.'
4918                   CALL message( 'define_netcdf_header', 'PA0251',             &
4919                                                                 0, 1, 0, 6, 0 )
4920                   extend = .FALSE.
4921                   RETURN
4922                ENDIF
4923             ELSE
4924                IF ( -1.0_wp /= netcdf_data(i) )  THEN
4925                   message_string = 'netCDF file for cross-sections ' //       &
4926                              TRIM( var ) // ' from previous run found,' //    &
4927                              ' but this file cannot be extended' //           &
4928                              ' due to mismatch in cross' //                   &
4929                              ' section levels.' //                            &
4930                              ' New file is created instead.'
4931                   CALL message( 'define_netcdf_header', 'PA0251',             &
4932                                                                 0, 1, 0, 6, 0 )
4933                   extend = .FALSE.
4934                   RETURN
4935                ENDIF
4936             ENDIF
4937          ENDDO
4938
4939          DEALLOCATE( netcdf_data )
4940
4941!
4942!--       Get the id of the time coordinate (unlimited coordinate) and its
4943!--       last index on the file. The next time level is pl2d..count+1.
4944!--       The current time must be larger than the last output time
4945!--       on the file.
4946          nc_stat = NF90_INQ_VARID( id_set_yz(av), 'time', id_var_time_yz(av) )
4947          CALL netcdf_handle_error( 'netcdf_define_header', 212 )
4948
4949          nc_stat = NF90_INQUIRE_VARIABLE( id_set_yz(av), id_var_time_yz(av), &
4950                                           dimids = id_dim_time_old )
4951          CALL netcdf_handle_error( 'netcdf_define_header', 213 )
4952          id_dim_time_yz(av) = id_dim_time_old(1)
4953
4954          nc_stat = NF90_INQUIRE_DIMENSION( id_set_yz(av), id_dim_time_yz(av), &
4955                                            len = ntime_count )
4956          CALL netcdf_handle_error( 'netcdf_define_header', 214 )
4957
4958!
4959!--       For non-parallel output use the last output time level of the netcdf
4960!--       file because the time dimension is unlimited. In case of parallel
4961!--       output the variable ntime_count could get the value of 9*10E36 because
4962!--       the time dimension is limited.
4963          IF ( netcdf_data_format < 5 ) do2d_yz_time_count(av) = ntime_count
4964
4965          nc_stat = NF90_GET_VAR( id_set_yz(av), id_var_time_yz(av),           &
4966                                  last_time_coordinate,                        &
4967                                  start = (/ do2d_yz_time_count(av) /),        &
4968                                  count = (/ 1 /) )
4969          CALL netcdf_handle_error( 'netcdf_define_header', 215 )
4970
4971          IF ( last_time_coordinate(1) >= simulated_time )  THEN
4972             message_string = 'netCDF file for cross sections ' //             &
4973                              TRIM( var ) // ' from previous run found,' //    &
4974                              '&but this file cannot be extended becaus' //    &
4975                              'e the current output time' //                   &
4976                              '&is less or equal than the last output t' //    &
4977                              'ime on this file.' //                           &
4978                              '&New file is created instead.'
4979             CALL message( 'define_netcdf_header', 'PA0252', 0, 1, 0, 6, 0 )
4980             do2d_yz_time_count(av) = 0
4981             extend = .FALSE.
4982             RETURN
4983          ENDIF
4984
4985          IF ( netcdf_data_format > 4 )