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

Last change on this file since 3435 was 3435, checked in by gronemeier, 3 years ago

new: terrain-following masked output; bugfixes: increase vertical dimension of gamma_w_green_sat by 1, add checks for masked output for chemistry_model and radiation_model, reordered calls to xxx_define_netcdf_grid in masked output part

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