Home

Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

netcdf_interface_mod.f90

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