Home

Context Navigation

netcdf_interface.f90 @ 1786

Last change on this file since 1786 was 1786, checked in by raasch, 7 years ago
pmc-change in server-client get-put, spectra-directives removed, spectra-package modularized
Property svn:keywords set to `Id`
File size: 218.8 KB

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