Home

Context Navigation

netcdf_interface_mod.f90 @ 3464

Last change on this file since 3464 was 3464, checked in by kanani, 6 years ago
from branch resler@3462: add MRT shaping function (radiation_model_mod), use basic constants (biometeorology_mod), adjust precision to wp (biometeorology_pt_mod)
Property svn:keywords set to `Id`
File size: 334.6 KB

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