Home

Context Navigation

netcdf_interface_mod.f90 @ 3187

Last change on this file since 3187 was 3187, checked in by sward, 7 years ago
Reworked agent pathfinding and modified output
Property svn:keywords set to `Id`
File size: 254.7 KB

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