Home

Context Navigation

netcdf_interface_mod.f90 @ 3525

Last change on this file since 3525 was 3525, checked in by kanani, 6 years ago
Changes related to clean-up of biometeorology (by dom_dwd_user)
Property svn:keywords set to `Id`
File size: 337.9 KB

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