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

Last change on this file since 3448 was 3448, checked in by kanani, 3 years ago

Implementation of human thermal indices (from branch biomet_p2 at r3444)

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