Home

Context Navigation

netcdf_interface_mod.f90 @ 3435

Last change on this file since 3435 was 3435, checked in by gronemeier, 6 years ago
new: terrain-following masked output; bugfixes: increase vertical dimension of gamma_w_green_sat by 1, add checks for masked output for chemistry_model and radiation_model, reordered calls to xxx_define_netcdf_grid in masked output part
Property svn:keywords set to `Id`
File size: 303.6 KB

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