Home

Context Navigation

netcdf_interface_mod.f90 @ 3467

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