Home

Context Navigation

netcdf_interface_mod.f90 @ 3459

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