Home

Context Navigation

netcdf_interface_mod.f90 @ 2012

Last change on this file since 2012 was 2012, checked in by kanani, 6 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 233.4 KB

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