Home

Context Navigation

netcdf_interface_mod.f90 @ 4227

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