Home

Context Navigation

netcdf_interface_mod.f90 @ 4196

Last change on this file since 4196 was 4196, checked in by gronemeier, 4 years ago

Consider rotation of model domain for claculation of Coriolis force:

check_parameters: Overwrite rotation_angle from namelist by value from static driver;
coriolis: Consider rotation of model domain;
header: Write information about rotation angle;
modules: Added rotation_angle;
netcdf_interface_mod: replaced rotation angle from input-netCDF file by namelist parameter 'rotation_angle';
ocean_mod: Consider rotation of model domain for calculating the Stokes drift;
parin: added rotation_angle to initialization_parameters;

Property svn:keywords set to Id

File size: 319.8 KB

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