Home

Context Navigation

netcdf_interface_mod.f90 @ 2239

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