Home

Context Navigation

source: palm/trunk/SOURCE/netcdf_interface_mod.f90 @ 2747

Last change on this file since 2747 was 2746, checked in by suehring, 7 years ago
Read information from statitic driver for resolved vegetation independently from land- or urban-surface model
Property svn:keywords set to `Id`
File size: 246.6 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 2746 2018-01-15 12:06:04Z suehring $
27	! Move flag plant canopy to modules
28	!
29	! 2718 2018-01-02 08:49:38Z maronga
30	! Corrected "Former revisions" section
31	!
32	! 2696 2017-12-14 17:12:51Z kanani
33	! Change in file header (GPL part)
34	! Implementation of uv exposure model (FK)
35	! Implemented checks for turbulence_closure_mod (TG)
36	! Implementation of chemistry module (FK)
37	! Bugfix in setting netcdf grids for LSM variables
38	! Enable setting of _FillValue attribute in output files (MS)
39	!
40	! 2512 2017-10-04 08:26:59Z raasch
41	! upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny
42	! no output of ghost layer data any more
43	!
44	! 2302 2017-07-03 14:07:20Z suehring
45	! Reading of 3D topography using NetCDF data type NC_BYTE
46	!
47	! 2292 2017-06-20 09:51:42Z schwenkel
48	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
49	! includes two more prognostic equations for cloud drop concentration (nc)
50	! and cloud water content (qc).
51	!
52	! 2270 2017-06-09 12:18:47Z maronga
53	! Removed 2 timeseries (shf_eb + qsws_eb). Removed _eb suffixes
54	!
55	! 2265 2017-06-08 16:58:28Z schwenkel
56	! Unused variables removed.
57	!
58	! 2239 2017-06-01 12:04:51Z suehring
59	! Bugfix xy-output of land-surface variables
60	!
61	! 2233 2017-05-30 18:08:54Z suehring
62	!
63	! 2232 2017-05-30 17:47:52Z suehring
64	! Adjustments to new topography and surface concept
65	!
66	! Topograpyh height arrays (zu_s_inner, zw_w_inner) are defined locally, output
67	! only if parallel netcdf.
68	!
69	! Build interface for topography input:
70	! - open file in read-only mode
71	! - read global attributes
72	! - read variables
73	!
74	! Bugfix in xy output (land-surface case)
75	!
76	! 2209 2017-04-19 09:34:46Z kanani
77	! Added support for plant canopy model output
78	!
79	! 2189 2017-03-21 09:29:52Z raasch
80	! bugfix: rho renamed rho_ocean for the cross section output
81	!
82	! 2109 2017-01-10 12:18:08Z raasch
83	! bugfix: length of character string netcdf_var_name extended to avoid problems
84	! which appeared in restart runs due to truncation
85	!
86	! 2040 2016-10-26 16:58:09Z gronemeier
87	! Increased number of possible statistic_regions to 99
88	!
89	! 2037 2016-10-26 11:15:40Z knoop
90	! Anelastic approximation implemented
91	!
92	! 2031 2016-10-21 15:11:58Z knoop
93	! renamed variable rho to rho_ocean
94	!
95	! 2011 2016-09-19 17:29:57Z kanani
96	! Flag urban_surface is now defined in module control_parameters,
97	! changed prefix for urban surface model output to "usm_",
98	! introduced control parameter varnamelength for LEN of trimvar.
99	!
100	! 2007 2016-08-24 15:47:17Z kanani
101	! Added support for new urban surface model (temporary modifications of
102	! SELECT CASE ( ) necessary, see variable trimvar),
103	! increased DIMENSION of do2d_unit, do3d_unit, id_var_do2d, id_var_do3d,
104	! increased LEN of char_cross_profiles, var_list, var_list_old
105	!
106	! 2000 2016-08-20 18:09:15Z knoop
107	! Forced header and separation lines into 80 columns
108	!
109	! 1990 2016-08-12 09:54:36Z gronemeier
110	! Bugfix: variable list was not written for time series output
111	!
112	! 1980 2016-07-29 15:51:57Z suehring
113	! Bugfix, in order to steer user-defined output, setting flag found explicitly
114	! to .F.
115	!
116	! 1976 2016-07-27 13:28:04Z maronga
117	! Removed remaining 2D land surface quantities. Definition of radiation
118	! quantities is now done directly in the respective module
119	!
120	! 1972 2016-07-26 07:52:02Z maronga
121	! Bugfix: wrong units for lsm quantities.
122	! Definition of grids for land surface quantities is now done directly in the
123	! respective module.
124	!
125	! 1960 2016-07-12 16:34:24Z suehring
126	! Additional labels and units for timeseries output of passive scalar-related
127	! quantities
128	!
129	! 1957 2016-07-07 10:43:48Z suehring
130	! flight module added
131	!
132	! 1850 2016-04-08 13:29:27Z maronga
133	! Module renamed
134	!
135	!
136	! 1833 2016-04-07 14:23:03Z raasch
137	! spectrum renamed spectra_mod
138	!
139	! 1786 2016-03-08 05:49:27Z raasch
140	! Bugfix: id_var_time_sp made public
141	!
142	! 1783 2016-03-06 18:36:17Z raasch
143	! netcdf interface has been modularized, former file netcdf renamed to
144	! netcdf_interface, creation of netcdf-dimensions and -variables moved to
145	! specific new subroutines create_netcdf_dim and create_netcdf_var,
146	! compression (deflation) of variables implemented,
147	! ibmy special cpp directive removed
148	!
149	! 1745 2016-02-05 13:06:51Z gronemeier
150	! Bugfix: recalculating ntdim_3d, ntdim_2d_xy/xz/yz when checking the
151	! extensibility of an existing file (only when using parallel NetCDF).
152	!
153	! 1691 2015-10-26 16:17:44Z maronga
154	! Added output of radiative heating rates for RRTMG. Corrected output of
155	! radiative fluxes
156	!
157	! 1682 2015-10-07 23:56:08Z knoop
158	! Code annotations made doxygen readable
159	!
160	! 1596 2015-05-21 09:34:28Z gronemeier
161	! Bugfix in masked data output. Read 'zu_3d' when trying to extend masked data
162	!
163	! 1551 2015-03-03 14:18:16Z maronga
164	! Added support for land surface model and radiation model output. In the course
165	! of this action a new vertical grid zs (soil) was introduced.
166	!
167	! 1353 2014-04-08 15:21:23Z heinze
168	! REAL constants provided with KIND-attribute
169	!
170	! 1322 2014-03-20 16:38:49Z raasch
171	! Forgotten ONLY-attribute added to USE-statements
172	!
173	! 1320 2014-03-20 08:40:49Z raasch
174	! ONLY-attribute added to USE-statements,
175	! kind-parameters added to all INTEGER and REAL declaration statements,
176	! kinds are defined in new module kinds,
177	! revision history before 2012 removed,
178	! comment fields (!:) to be used for variable explanations added to
179	! all variable declaration statements
180	!
181	! 1308 2014-03-13 14:58:42Z fricke
182	! +ntime_count, oldmode
183	! Adjust NF90_CREATE and NF90_OPEN statement for parallel output
184	! To increase the performance for parallel output, the following is done:
185	! - Limit time dimension
186	! - Values of axis data are only written by PE0
187	! - No fill is set for all variables
188	! Check the number of output time levels for restart jobs
189	!
190	! 1206 2013-07-18 12:49:16Z witha
191	! Bugfix: typo in preprocessor directive in subroutine open_write_netcdf_file
192	!
193	! 1092 2013-02-02 11:24:22Z raasch
194	! unused variables removed
195	!
196	! 1053 2012-11-13 17:11:03Z hoffmann
197	! +qr, nr, prr
198	!
199	! 1036 2012-10-22 13:43:42Z raasch
200	! code put under GPL (PALM 3.9)
201	!
202	! 1031 2012-10-19 14:35:30Z raasch
203	! netCDF4 without parallel file support implemented, new routines
204	! create_netcdf_file and open_write_netcdf_file at end
205	!
206	! 992 2012-09-05 15:08:26Z hoffmann
207	! Removal of the informative messages PA0352 and PA0353.
208	!
209	! 983 2012-08-21 14:17:57Z hoffmann
210	! Bugfix in cross_profiles.
211	!
212	! 964 2012-07-26 09:14:24Z raasch
213	! rev 951 and 959 reformatted
214	!
215	! 959 2012-07-24 13:13:41Z hoffmann
216	! Bugfix in cross_profiles. It is not allowed to arrange more than 100
217	! profiles with cross_profiles.
218	!
219	! 951 2012-07-19 14:22:52Z hoffmann
220	! cross_profiles, profile_rows, profile_columns are written to netCDF header
221	!
222	! Revision 1.1 2005/05/18 15:37:16 raasch
223	! Initial revision
224	!
225	!
226	! Description:
227	! ------------
228	!> In case of extend = .FALSE.:
229	!> Define all necessary dimensions, axes and variables for the different
230	!> netCDF datasets. This subroutine is called from check_open after a new
231	!> dataset is created. It leaves the open netCDF files ready to write.
232	!>
233	!> In case of extend = .TRUE.:
234	!> Find out if dimensions and variables of an existing file match the values
235	!> of the actual run. If so, get all necessary information (ids, etc.) from
236	!> this file.
237	!>
238	!> Parameter av can assume values 0 (non-averaged data) and 1 (time averaged
239	!> data)
240	!>
241	!> @todo calculation of output time levels for parallel NetCDF still does not
242	!> cover every exception (change of dt_do, end_time in restart)
243	!> @todo timeseries and profile output still needs to be rewritten to allow
244	!> modularization
245	!------------------------------------------------------------------------------!
246	MODULE netcdf_interface
247
248	USE control_parameters, &
249	ONLY: max_masks, fl_max, var_fl_max, varnamelength
250	USE kinds
251	#if defined( __netcdf )
252	USE NETCDF
253	#endif
254
255	PRIVATE
256
257	INTEGER(iwp), PARAMETER :: dopr_norm_num = 7, dopts_num = 29, dots_max = 100
258
259	CHARACTER (LEN=6), DIMENSION(dopr_norm_num) :: dopr_norm_names = &
260	(/ 'wpt0 ', 'ws2 ', 'tsw2 ', 'ws3 ', 'ws2tsw', 'wstsw2', &
261	'z_i ' /)
262
263	CHARACTER (LEN=6), DIMENSION(dopr_norm_num) :: dopr_norm_longnames = &
264	(/ 'wpt0 ', 'w2 ', 'tw2 ', 'w3 ', 'w2tw', 'wt*w2', &
265	'z_i ' /)
266
267	CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_label = &
268	(/ 'tnpt ', 'x_ ', 'y_ ', 'z_ ', 'z_abs ', 'u ', &
269	'v ', 'w ', 'u" ', 'v" ', 'w" ', 'npt_up ', &
270	'w_up ', 'w_down ', 'radius ', 'r_min ', 'r_max ', 'npt_max', &
271	'npt_min', 'x2 ', 'y2 ', 'z2 ', 'u2 ', 'v*2 ', &
272	'w2 ', 'u"2 ', 'v"2 ', 'w"2 ', 'npt2 ' /)
273
274	CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_unit = &
275	(/ 'number ', 'm ', 'm ', 'm ', 'm ', 'm/s ', &
276	'm/s ', 'm/s ', 'm/s ', 'm/s ', 'm/s ', 'number ', &
277	'm/s ', 'm/s ', 'm ', 'm ', 'm ', 'number ', &
278	'number ', 'm2 ', 'm2 ', 'm2 ', 'm2/s2 ', 'm2/s2 ', &
279	'm2/s2 ', 'm2/s2 ', 'm2/s2 ', 'm2/s2 ', 'number2' /)
280
281	INTEGER(iwp) :: dots_num = 25 !< number of timeseries defined by default
282	INTEGER(iwp) :: dots_soil = 26 !< starting index for soil-timeseries
283	INTEGER(iwp) :: dots_rad = 32 !< starting index for radiation-timeseries
284
285	CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label = &
286	(/ 'E ', 'E* ', 'dt ', &
287	'u* ', 'th* ', 'umax ', &
288	'vmax ', 'wmax ', 'div_new ', &
289	'div_old ', 'z_i_wpt ', 'z_i_pt ', &
290	'w* ', 'w"pt"0 ', 'w"pt" ', &
291	'wpt ', 'pt(0) ', 'pt(z_mo) ', &
292	'w"u"0 ', 'w"v"0 ', 'w"q"0 ', &
293	'ol ', 'q* ', 'w"s" ', &
294	's* ', 'ghf ', 'qsws_liq ', &
295	'qsws_soil ', 'qsws_veg ', 'r_a ', &
296	'r_s ', &
297	'rad_net ', 'rad_lw_in ', 'rad_lw_out ', &
298	'rad_sw_in ', 'rad_sw_out ', 'rrtm_aldif ', &
299	'rrtm_aldir ', 'rrtm_asdif ', 'rrtm_asdir ', &
300	( 'unknown ', i9 = 1, dots_max-40 ) /)
301
302	CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit = &
303	(/ 'm2/s2 ', 'm2/s2 ', 's ', &
304	'm/s ', 'K ', 'm/s ', &
305	'm/s ', 'm/s ', 's-1 ', &
306	's-1 ', 'm ', 'm ', &
307	'm/s ', 'K m/s ', 'K m/s ', &
308	'K m/s ', 'K ', 'K ', &
309	'm2/s2 ', 'm2/s2 ', 'kg m/s ', &
310	'm ', 'kg/kg ', 'kg m/(kg s) ', &
311	'kg/kg ', 'W/m2 ', 'W/m2 ', &
312	'W/m2 ', 'W/m2 ', 's/m ', &
313	's/m ', &
314	'W/m2 ', 'W/m2 ', 'W/m2 ', &
315	'W/m2 ', 'W/m2 ', ' ', &
316	' ', ' ', ' ', &
317	( 'unknown ', i9 = 1, dots_max-40 ) /)
318
319	CHARACTER (LEN=16) :: heatflux_output_unit !< unit for heatflux output
320	CHARACTER (LEN=16) :: waterflux_output_unit !< unit for waterflux output
321	CHARACTER (LEN=16) :: momentumflux_output_unit !< unit for momentumflux output
322
323	CHARACTER (LEN=9), DIMENSION(300) :: dopr_unit = 'unknown'
324
325	CHARACTER (LEN=7), DIMENSION(0:1,500) :: do2d_unit, do3d_unit
326
327	CHARACTER (LEN=16), DIMENSION(25) :: prt_var_names = &
328	(/ 'pt_age ', 'pt_dvrp_size ', 'pt_origin_x ', &
329	'pt_origin_y ', 'pt_origin_z ', 'pt_radius ', &
330	'pt_speed_x ', 'pt_speed_y ', 'pt_speed_z ', &
331	'pt_weight_factor', 'pt_x ', 'pt_y ', &
332	'pt_z ', 'pt_color ', 'pt_group ', &
333	'pt_tailpoints ', 'pt_tail_id ', 'pt_density_ratio', &
334	'pt_exp_arg ', 'pt_exp_term ', 'not_used ', &
335	'not_used ', 'not_used ', 'not_used ', &
336	'not_used ' /)
337
338	CHARACTER (LEN=16), DIMENSION(25) :: prt_var_units = &
339	(/ 'seconds ', 'meters ', 'meters ', &
340	'meters ', 'meters ', 'meters ', &
341	'm/s ', 'm/s ', 'm/s ', &
342	'factor ', 'meters ', 'meters ', &
343	'meters ', 'none ', 'none ', &
344	'none ', 'none ', 'ratio ', &
345	'none ', 'none ', 'not_used ', &
346	'not_used ', 'not_used ', 'not_used ', &
347	'not_used ' /)
348
349	CHARACTER(LEN=20), DIMENSION(11) :: netcdf_precision = ' '
350	CHARACTER(LEN=40) :: netcdf_data_format_string
351
352	INTEGER(iwp) :: id_dim_prtnum, id_dim_time_fl, id_dim_time_pr, id_dim_time_prt, &
353	id_dim_time_pts, id_dim_time_sp, id_dim_time_ts, id_dim_x_sp, &
354	id_dim_y_sp, id_dim_zu_sp, id_dim_zw_sp, id_set_fl, id_set_pr, &
355	id_set_prt, id_set_pts, id_set_sp, id_set_ts, id_var_time_fl, &
356	id_var_prtnum, id_var_rnop_prt, id_var_time_pr, id_var_time_prt, &
357	id_var_time_pts, id_var_time_sp, id_var_time_ts, id_var_x_sp, &
358	id_var_y_sp, id_var_zu_sp, id_var_zw_sp, nc_stat
359
360
361	INTEGER(iwp), DIMENSION(0:1) :: id_dim_time_xy, id_dim_time_xz, &
362	id_dim_time_yz, id_dim_time_3d, id_dim_x_xy, id_dim_xu_xy, &
363	id_dim_x_xz, id_dim_xu_xz, id_dim_x_yz, id_dim_xu_yz, &
364	id_dim_x_3d, id_dim_xu_3d, id_dim_y_xy, id_dim_yv_xy, &
365	id_dim_y_xz, id_dim_yv_xz, id_dim_y_yz, id_dim_yv_yz, &
366	id_dim_y_3d, id_dim_yv_3d, id_dim_zs_xy, id_dim_zs_xz, &
367	id_dim_zs_yz, id_dim_zs_3d, id_dim_zu_xy, id_dim_zu1_xy, &
368	id_dim_zu_xz, id_dim_zu_yz, id_dim_zu_3d, id_dim_zw_xy, &
369	id_dim_zw_xz, id_dim_zw_yz, id_dim_zw_3d, id_set_xy, &
370	id_set_xz, id_set_yz, id_set_3d, id_var_ind_x_yz, &
371	id_var_ind_y_xz, id_var_ind_z_xy, id_var_time_xy, &
372	id_var_time_xz, id_var_time_yz, id_var_time_3d, id_var_x_xy, &
373	id_var_xu_xy, id_var_x_xz, id_var_xu_xz, id_var_x_yz, &
374	id_var_xu_yz, id_var_x_3d, id_var_xu_3d, id_var_y_xy, &
375	id_var_yv_xy, id_var_y_xz, id_var_yv_xz, id_var_y_yz, &
376	id_var_yv_yz, id_var_y_3d, id_var_yv_3d, id_var_zs_xy, &
377	id_var_zs_xz, id_var_zs_yz, id_var_zs_3d, id_var_zusi_xy, &
378	id_var_zusi_3d, id_var_zu_xy, id_var_zu1_xy, id_var_zu_xz, &
379	id_var_zu_yz, id_var_zu_3d, id_var_zwwi_xy, id_var_zwwi_3d, &
380	id_var_zw_xy, id_var_zw_xz, id_var_zw_yz, id_var_zw_3d
381
382	INTEGER :: netcdf_data_format = 2 !< NetCDF3 64bit offset format
383	INTEGER :: netcdf_deflate = 0 !< NetCDF compression, default: no
384	!< compression
385
386	INTEGER(iwp) :: dofl_time_count
387	INTEGER(iwp), DIMENSION(10) :: id_var_dospx, id_var_dospy
388	INTEGER(iwp), DIMENSION(20) :: id_var_prt
389	INTEGER(iwp), DIMENSION(11) :: nc_precision
390	INTEGER(iwp), DIMENSION(dopr_norm_num) :: id_var_norm_dopr
391
392	INTEGER(iwp), DIMENSION(fl_max) :: id_dim_x_fl, id_dim_y_fl, id_dim_z_fl
393	INTEGER(iwp), DIMENSION(fl_max) :: id_var_x_fl, id_var_y_fl, id_var_z_fl
394
395	CHARACTER (LEN=20), DIMENSION(fl_max*var_fl_max) :: dofl_label
396	CHARACTER (LEN=20), DIMENSION(fl_max*var_fl_max) :: dofl_unit
397	CHARACTER (LEN=20), DIMENSION(fl_max) :: dofl_dim_label_x
398	CHARACTER (LEN=20), DIMENSION(fl_max) :: dofl_dim_label_y
399	CHARACTER (LEN=20), DIMENSION(fl_max) :: dofl_dim_label_z
400
401	INTEGER(iwp), DIMENSION(fl_max*var_fl_max) :: id_var_dofl
402
403	INTEGER(iwp), DIMENSION(dopts_num,0:10) :: id_var_dopts
404	INTEGER(iwp), DIMENSION(0:1,500) :: id_var_do2d, id_var_do3d
405	INTEGER(iwp), DIMENSION(100,0:99) :: id_dim_z_pr, id_var_dopr, &
406	id_var_z_pr
407	INTEGER(iwp), DIMENSION(dots_max,0:99) :: id_var_dots
408
409	!
410	!-- Masked output
411	CHARACTER (LEN=7), DIMENSION(max_masks,0:1,100) :: domask_unit
412
413	LOGICAL :: output_for_t0 = .FALSE.
414
415	INTEGER(iwp), DIMENSION(1:max_masks,0:1) :: id_dim_time_mask, id_dim_x_mask, &
416	id_dim_xu_mask, id_dim_y_mask, id_dim_yv_mask, id_dim_zs_mask, &
417	id_dim_zu_mask, id_dim_zw_mask, &
418	id_set_mask, &
419	id_var_time_mask, id_var_x_mask, id_var_xu_mask, &
420	id_var_y_mask, id_var_yv_mask, id_var_zs_mask, &
421	id_var_zu_mask, id_var_zw_mask, &
422	id_var_zusi_mask, id_var_zwwi_mask
423
424	INTEGER(iwp), DIMENSION(1:max_masks,0:1,100) :: id_var_domask
425
426	REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute
427
428
429	PUBLIC dofl_dim_label_x, dofl_dim_label_y, dofl_dim_label_z, dofl_label, &
430	dofl_time_count, dofl_unit, domask_unit, dopr_unit, dopts_num, &
431	dots_label, dots_max, dots_num, dots_rad, dots_soil, dots_unit, &
432	do2d_unit, do3d_unit, fill_value, &
433	id_set_fl, id_set_mask, id_set_pr, &
434	id_set_prt, id_set_pts, id_set_sp, id_set_ts, id_set_xy, id_set_xz,&
435	id_set_yz, id_set_3d, id_var_domask, id_var_dofl, id_var_dopr, &
436	id_var_dopts, id_var_dospx, id_var_dospy, id_var_dots, id_var_do2d,&
437	id_var_do3d, id_var_norm_dopr, id_var_time_fl, id_var_time_mask, &
438	id_var_time_pr, id_var_time_pts, id_var_time_sp, id_var_time_ts, &
439	id_var_time_xy, id_var_time_xz, id_var_time_yz, id_var_time_3d, &
440	id_var_x_fl, id_var_y_fl, id_var_z_fl, nc_stat, &
441	netcdf_data_format, netcdf_data_format_string, netcdf_deflate, &
442	netcdf_precision, output_for_t0, heatflux_output_unit, &
443	waterflux_output_unit, momentumflux_output_unit
444
445	SAVE
446
447	INTERFACE netcdf_create_dim
448	MODULE PROCEDURE netcdf_create_dim
449	END INTERFACE netcdf_create_dim
450
451	INTERFACE netcdf_create_file
452	MODULE PROCEDURE netcdf_create_file
453	END INTERFACE netcdf_create_file
454
455	INTERFACE netcdf_create_var
456	MODULE PROCEDURE netcdf_create_var
457	END INTERFACE netcdf_create_var
458
459	INTERFACE netcdf_define_header
460	MODULE PROCEDURE netcdf_define_header
461	END INTERFACE netcdf_define_header
462
463	INTERFACE netcdf_handle_error
464	MODULE PROCEDURE netcdf_handle_error
465	END INTERFACE netcdf_handle_error
466
467	INTERFACE netcdf_open_write_file
468	MODULE PROCEDURE netcdf_open_write_file
469	END INTERFACE netcdf_open_write_file
470
471	PUBLIC netcdf_create_file, netcdf_define_header, &
472	netcdf_handle_error, netcdf_open_write_file
473
474	CONTAINS
475
476	SUBROUTINE netcdf_define_header( callmode, extend, av )
477
478	#if defined( __netcdf )
479
480	USE arrays_3d, &
481	ONLY: zu, zw
482
483	#if defined( __chem )
484	USE chemistry_model_mod, &
485	ONLY: chem_define_netcdf_grid
486	#endif
487
488	USE constants, &
489	ONLY: pi
490
491	USE control_parameters, &
492	ONLY: air_chemistry, averaging_interval, averaging_interval_pr, &
493	data_output_pr, domask, dopr_n, &
494	dopr_time_count, dopts_time_count, dots_time_count, &
495	do2d, do2d_xz_time_count, do3d, &
496	do2d_yz_time_count, dt_data_output_av, dt_do2d_xy, dt_do2d_xz, &
497	dt_do2d_yz, dt_do3d, mask_size, do2d_xy_time_count, &
498	do3d_time_count, domask_time_count, end_time, land_surface, &
499	mask_size_l, mask_i, mask_i_global, mask_j, mask_j_global, &
500	mask_k_global, message_string, mid, ntdim_2d_xy, &
501	ntdim_2d_xz, ntdim_2d_yz, ntdim_3d, nz_do3d, plant_canopy, &
502	prt_time_count, run_description_header, section, simulated_time,&
503	simulated_time_at_begin, skip_time_data_output_av, &
504	skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz, &
505	skip_time_do3d, topography, num_leg, num_var_fl, &
506	urban_surface, uv_exposure
507
508	USE grid_variables, &
509	ONLY: dx, dy, zu_s_inner, zw_w_inner
510
511	USE indices, &
512	ONLY: nx, nxl, nxr, ny, nys, nyn, nz ,nzb, nzt
513
514	USE kinds
515
516	USE land_surface_model_mod, &
517	ONLY: lsm_define_netcdf_grid, nzb_soil, nzt_soil, nzs, zs
518
519	USE pegrid
520
521	USE particle_attributes, &
522	ONLY: number_of_particle_groups
523
524	USE plant_canopy_model_mod, &
525	ONLY: pcm_define_netcdf_grid
526
527	USE profil_parameter, &
528	ONLY: crmax, cross_profiles, dopr_index, profile_columns, profile_rows
529
530	USE radiation_model_mod, &
531	ONLY: radiation, radiation_define_netcdf_grid
532
533	USE spectra_mod, &
534	ONLY: averaging_interval_sp, comp_spectra_level, data_output_sp, dosp_time_count, spectra_direction
535
536	USE statistics, &
537	ONLY: hom, statistic_regions
538
539	USE turbulence_closure_mod, &
540	ONLY: tcm_define_netcdf_grid
541
542	USE urban_surface_mod, &
543	ONLY: usm_define_netcdf_grid
544
545	USE uv_exposure_model_mod, &
546	ONLY: uvem_define_netcdf_grid
547
548
549	IMPLICIT NONE
550
551	CHARACTER (LEN=3) :: suffix !<
552	CHARACTER (LEN=2), INTENT (IN) :: callmode !<
553	CHARACTER (LEN=4) :: grid_x !<
554	CHARACTER (LEN=4) :: grid_y !<
555	CHARACTER (LEN=4) :: grid_z !<
556	CHARACTER (LEN=6) :: mode !<
557	CHARACTER (LEN=10) :: precision !<
558	CHARACTER (LEN=10) :: var !<
559	CHARACTER (LEN=20) :: netcdf_var_name !<
560	CHARACTER (LEN=varnamelength) :: trimvar !< TRIM of output-variable string
561	CHARACTER (LEN=80) :: time_average_text !<
562	CHARACTER (LEN=4000) :: char_cross_profiles !<
563	CHARACTER (LEN=4000) :: var_list !<
564	CHARACTER (LEN=4000) :: var_list_old !<
565
566	CHARACTER (LEN=100), DIMENSION(1:crmax) :: cross_profiles_adj !<
567	CHARACTER (LEN=100), DIMENSION(1:crmax) :: cross_profiles_char !<
568
569	INTEGER(iwp) :: av !<
570	INTEGER(iwp) :: cross_profiles_count !<
571	INTEGER(iwp) :: cross_profiles_maxi !<
572	INTEGER(iwp) :: delim !<
573	INTEGER(iwp) :: delim_old !<
574	INTEGER(iwp) :: file_id !<
575	INTEGER(iwp) :: i !<
576	INTEGER(iwp) :: id_last !<
577	INTEGER(iwp) :: id_x !<
578	INTEGER(iwp) :: id_y !<
579	INTEGER(iwp) :: id_z !<
580	INTEGER(iwp) :: j !<
581	INTEGER(iwp) :: k !<
582	INTEGER(iwp) :: kk !<
583	INTEGER(iwp) :: ns !<
584	INTEGER(iwp) :: ns_do !< actual value of ns for soil model data
585	INTEGER(iwp) :: ns_old !<
586	INTEGER(iwp) :: ntime_count !< number of time levels found in file
587	INTEGER(iwp) :: nz_old !<
588	INTEGER(iwp) :: l !<
589
590	INTEGER(iwp), SAVE :: oldmode !<
591
592	INTEGER(iwp), DIMENSION(1) :: id_dim_time_old !<
593	INTEGER(iwp), DIMENSION(1) :: id_dim_x_yz_old !<
594	INTEGER(iwp), DIMENSION(1) :: id_dim_y_xz_old !<
595	INTEGER(iwp), DIMENSION(1) :: id_dim_zu_sp_old !<
596	INTEGER(iwp), DIMENSION(1) :: id_dim_zu_xy_old !<
597	INTEGER(iwp), DIMENSION(1) :: id_dim_zu_3d_old !<
598	INTEGER(iwp), DIMENSION(1) :: id_dim_zu_mask_old !<
599
600
601	INTEGER(iwp), DIMENSION(1:crmax) :: cross_profiles_numb !<
602
603	LOGICAL :: found !<
604
605	LOGICAL, INTENT (INOUT) :: extend !<
606
607	LOGICAL, SAVE :: init_netcdf = .FALSE. !<
608
609	REAL(wp), DIMENSION(1) :: last_time_coordinate !< last time value in file
610
611	REAL(wp), DIMENSION(:), ALLOCATABLE :: netcdf_data !<
612	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: netcdf_data_2d !<
613
614
615	!
616	!-- Initializing actions
617	IF ( .NOT. init_netcdf ) THEN
618	!
619	!-- Check and set accuracy for netCDF output. First set default value
620	nc_precision = NF90_REAL4
621
622	i = 1
623	DO WHILE ( netcdf_precision(i) /= ' ' )
624	j = INDEX( netcdf_precision(i), '_' )
625	IF ( j == 0 ) THEN
626	WRITE ( message_string, * ) 'netcdf_precision must contain a ', &
627	'"_"netcdf_precision(', i, ')="', &
628	TRIM( netcdf_precision(i) ),'"'
629	CALL message( 'netcdf_define_header', 'PA0241', 2, 2, 0, 6, 0 )
630	ENDIF
631
632	var = netcdf_precision(i)(1:j-1)
633	precision = netcdf_precision(i)(j+1:)
634
635	IF ( precision == 'NF90_REAL4' ) THEN
636	j = NF90_REAL4
637	ELSEIF ( precision == 'NF90_REAL8' ) THEN
638	j = NF90_REAL8
639	ELSE
640	WRITE ( message_string, * ) 'illegal netcdf precision: ', &
641	'netcdf_precision(', i, ')="', &
642	TRIM( netcdf_precision(i) ),'"'
643	CALL message( 'netcdf_define_header', 'PA0242', 1, 2, 0, 6, 0 )
644	ENDIF
645
646	SELECT CASE ( var )
647	CASE ( 'xy' )
648	nc_precision(1) = j
649	CASE ( 'xz' )
650	nc_precision(2) = j
651	CASE ( 'yz' )
652	nc_precision(3) = j
653	CASE ( '2d' )
654	nc_precision(1:3) = j
655	CASE ( '3d' )
656	nc_precision(4) = j
657	CASE ( 'pr' )
658	nc_precision(5) = j
659	CASE ( 'ts' )
660	nc_precision(6) = j
661	CASE ( 'sp' )
662	nc_precision(7) = j
663	CASE ( 'prt' )
664	nc_precision(8) = j
665	CASE ( 'masks' )
666	nc_precision(11) = j
667	CASE ( 'fl' )
668	nc_precision(9) = j
669	CASE ( 'all' )
670	nc_precision = j
671
672	CASE DEFAULT
673	WRITE ( message_string, * ) 'unknown variable in inipar ', &
674	'assignment: netcdf_precision(', i, ')="', &
675	TRIM( netcdf_precision(i) ),'"'
676	CALL message( 'netcdf_define_header', 'PA0243', 1, 2, 0, 6, 0 )
677
678	END SELECT
679
680	i = i + 1
681	IF ( i > 50 ) EXIT
682	ENDDO
683
684	!
685	!-- Check for allowed parameter range
686	IF ( netcdf_deflate < 0 .OR. netcdf_deflate > 9 ) THEN
687	WRITE ( message_string, '(A,I3,A)' ) 'netcdf_deflate out of ' // &
688	'range & given value: ', netcdf_deflate, &
689	', allowed range: 0-9'
690	CALL message( 'netcdf_define_header', 'PA0355', 2, 2, 0, 6, 0 )
691	ENDIF
692	!
693	!-- Data compression does not work with parallel NetCDF/HDF5
694	IF ( netcdf_deflate > 0 .AND. netcdf_data_format /= 3 ) THEN
695	message_string = 'netcdf_deflate reset to 0'
696	CALL message( 'netcdf_define_header', 'PA0356', 0, 1, 0, 6, 0 )
697
698	netcdf_deflate = 0
699	ENDIF
700
701	init_netcdf = .TRUE.
702
703	ENDIF
704
705	!
706	!-- Determine the mode to be processed
707	IF ( extend ) THEN
708	mode = callmode // '_ext'
709	ELSE
710	mode = callmode // '_new'
711	ENDIF
712
713	!
714	!-- Select the mode to be processed. Possibilities are 3d, mask, xy, xz, yz,
715	!-- pr and ts.
716	SELECT CASE ( mode )
717
718	CASE ( 'ma_new' )
719
720	!
721	!-- decompose actual parameter file_id (=formal parameter av) into
722	!-- mid and av
723	file_id = av
724	IF ( file_id <= 200+max_masks ) THEN
725	mid = file_id - 200
726	av = 0
727	ELSE
728	mid = file_id - (200+max_masks)
729	av = 1
730	ENDIF
731
732	!
733	!-- Define some global attributes of the dataset
734	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
735	'Conventions', 'COARDS' )
736	CALL netcdf_handle_error( 'netcdf_define_header', 464 )
737
738	IF ( av == 0 ) THEN
739	time_average_text = ' '
740	ELSE
741	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
742	averaging_interval
743	ENDIF
744	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'title', &
745	TRIM( run_description_header ) // &
746	TRIM( time_average_text ) )
747	CALL netcdf_handle_error( 'netcdf_define_header', 465 )
748	IF ( av == 1 ) THEN
749	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval
750	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
751	'time_avg', TRIM( time_average_text ) )
752	CALL netcdf_handle_error( 'netcdf_define_header', 466 )
753	ENDIF
754
755	!
756	!-- Define time coordinate for volume data (unlimited dimension)
757	CALL netcdf_create_dim( id_set_mask(mid,av), 'time', NF90_UNLIMITED, &
758	id_dim_time_mask(mid,av), 467 )
759	CALL netcdf_create_var( id_set_mask(mid,av), &
760	(/ id_dim_time_mask(mid,av) /), 'time', &
761	NF90_DOUBLE, id_var_time_mask(mid,av), &
762	'seconds', '', 468, 469, 000 )
763	!
764	!-- Define spatial dimensions and coordinates:
765	!-- Define vertical coordinate grid (zu grid)
766	CALL netcdf_create_dim( id_set_mask(mid,av), 'zu_3d', &
767	mask_size(mid,3), id_dim_zu_mask(mid,av), &
768	470 )
769	CALL netcdf_create_var( id_set_mask(mid,av), &
770	(/ id_dim_zu_mask(mid,av) /), 'zu_3d', &
771	NF90_DOUBLE, id_var_zu_mask(mid,av), &
772	'meters', '', 471, 472, 000 )
773	!
774	!-- Define vertical coordinate grid (zw grid)
775	CALL netcdf_create_dim( id_set_mask(mid,av), 'zw_3d', &
776	mask_size(mid,3), id_dim_zw_mask(mid,av), &
777	473 )
778	CALL netcdf_create_var( id_set_mask(mid,av), &
779	(/ id_dim_zw_mask(mid,av) /), 'zw_3d', &
780	NF90_DOUBLE, id_var_zw_mask(mid,av), &
781	'meters', '', 474, 475, 000 )
782	!
783	!-- Define x-axis (for scalar position)
784	CALL netcdf_create_dim( id_set_mask(mid,av), 'x', mask_size(mid,1), &
785	id_dim_x_mask(mid,av), 476 )
786	CALL netcdf_create_var( id_set_mask(mid,av), &
787	(/ id_dim_x_mask(mid,av) /), 'x', &
788	NF90_DOUBLE, id_var_x_mask(mid,av), &
789	'meters', '', 477, 478, 000 )
790	!
791	!-- Define x-axis (for u position)
792	CALL netcdf_create_dim( id_set_mask(mid,av), 'xu', mask_size(mid,1), &
793	id_dim_xu_mask(mid,av), 479 )
794	CALL netcdf_create_var( id_set_mask(mid,av), &
795	(/ id_dim_xu_mask(mid,av) /), 'xu', &
796	NF90_DOUBLE, id_var_xu_mask(mid,av), &
797	'meters', '', 480, 481, 000 )
798	!
799	!-- Define y-axis (for scalar position)
800	CALL netcdf_create_dim( id_set_mask(mid,av), 'y', mask_size(mid,2), &
801	id_dim_y_mask(mid,av), 482 )
802	CALL netcdf_create_var( id_set_mask(mid,av), &
803	(/ id_dim_y_mask(mid,av) /), 'y', &
804	NF90_DOUBLE, id_var_y_mask(mid,av), &
805	'meters', '', 483, 484, 000 )
806	!
807	!-- Define y-axis (for v position)
808	CALL netcdf_create_dim( id_set_mask(mid,av), 'yv', mask_size(mid,2), &
809	id_dim_yv_mask(mid,av), 485 )
810	CALL netcdf_create_var( id_set_mask(mid,av), &
811	(/ id_dim_yv_mask(mid,av) /), &
812	'yv', NF90_DOUBLE, id_var_yv_mask(mid,av), &
813	'meters', '', 486, 487, 000 )
814	!
815	!-- In case of non-flat topography define 2d-arrays containing the height
816	!-- information. Only for parallel netcdf output.
817	IF ( TRIM( topography ) /= 'flat' .AND. &
818	netcdf_data_format > 4 ) THEN
819	!
820	!-- Define zusi = zu(nzb_s_inner)
821	CALL netcdf_create_var( id_set_mask(mid,av), &
822	(/ id_dim_x_mask(mid,av), &
823	id_dim_y_mask(mid,av) /), 'zusi', &
824	NF90_DOUBLE, id_var_zusi_mask(mid,av), &
825	'meters', 'zu(nzb_s_inner)', 488, 489, &
826	490 )
827	!
828	!-- Define zwwi = zw(nzb_w_inner)
829	CALL netcdf_create_var( id_set_mask(mid,av), &
830	(/ id_dim_x_mask(mid,av), &
831	id_dim_y_mask(mid,av) /), 'zwwi', &
832	NF90_DOUBLE, id_var_zwwi_mask(mid,av), &
833	'meters', 'zw(nzb_w_inner)', 491, 492, &
834	493 )
835	ENDIF
836
837	IF ( land_surface ) THEN
838	!
839	!-- Define vertical coordinate grid (zw grid)
840	CALL netcdf_create_dim( id_set_mask(mid,av), 'zs_3d', &
841	mask_size(mid,3), id_dim_zs_mask(mid,av), &
842	536 )
843	CALL netcdf_create_var( id_set_mask(mid,av), &
844	(/ id_dim_zs_mask(mid,av) /), 'zs_3d', &
845	NF90_DOUBLE, id_var_zs_mask(mid,av), &
846	'meters', '', 537, 555, 000 )
847	ENDIF
848
849	!
850	!-- Define the variables
851	var_list = ';'
852	i = 1
853
854	DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
855	!
856	!-- Temporary solution to account for data output within the new urban
857	!-- surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
858	trimvar = TRIM( domask(mid,av,i) )
859	IF ( urban_surface .AND. trimvar(1:4) == 'usm_' ) THEN
860	trimvar = 'usm_output'
861	ENDIF
862	!
863	!-- Check for the grid
864	found = .FALSE.
865	SELECT CASE ( trimvar )
866	!
867	!-- Most variables are defined on the scalar grid
868	CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', &
869	'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', &
870	'rho_ocean', 's', 'sa', 'vpt' )
871
872	grid_x = 'x'
873	grid_y = 'y'
874	grid_z = 'zu'
875	!
876	!-- u grid
877	CASE ( 'u' )
878
879	grid_x = 'xu'
880	grid_y = 'y'
881	grid_z = 'zu'
882	!
883	!-- v grid
884	CASE ( 'v' )
885
886	grid_x = 'x'
887	grid_y = 'yv'
888	grid_z = 'zu'
889	!
890	!-- w grid
891	CASE ( 'w' )
892
893	grid_x = 'x'
894	grid_y = 'y'
895	grid_z = 'zw'
896
897	!
898	!-- Block of urban surface model outputs
899	CASE ( 'usm_output' )
900
901	CALL usm_define_netcdf_grid( domask(mid,av,i), found, &
902	grid_x, grid_y, grid_z )
903
904	CASE DEFAULT
905
906	CALL tcm_define_netcdf_grid( domask(mid,av,i), found, &
907	grid_x, grid_y, grid_z )
908
909	!
910	!-- Check for land surface quantities
911	IF ( land_surface ) THEN
912	CALL lsm_define_netcdf_grid( domask(mid,av,i), found, &
913	grid_x, grid_y, grid_z )
914	ENDIF
915	!
916	!-- Check for plant canopy quantities
917	IF ( .NOT. found .AND. plant_canopy ) THEN
918	CALL pcm_define_netcdf_grid( domask(mid,av,i), found, &
919	grid_x, grid_y, grid_z )
920	ENDIF
921
922	!
923	!-- Check for radiation quantities
924	IF ( .NOT. found .AND. radiation ) THEN
925	CALL radiation_define_netcdf_grid( domask(mid,av,i), &
926	found, grid_x, grid_y,&
927	grid_z )
928	ENDIF
929
930	!
931	!-- Check for chemistry quantities
932	#if defined( __chem )
933	IF ( .NOT. found .AND. air_chemistry ) THEN
934	CALL chem_define_netcdf_grid( domask(mid,av,i), &
935	found, grid_x, grid_y, &
936	grid_z )
937	ENDIF
938	#endif
939
940	!
941	!-- Check for user-defined quantities
942	IF ( .NOT. found ) THEN
943	CALL user_define_netcdf_grid( domask(mid,av,i), found, &
944	grid_x, grid_y, grid_z )
945	ENDIF
946
947	IF ( .NOT. found ) THEN
948	WRITE ( message_string, * ) 'no grid defined for', &
949	' variable ', TRIM( domask(mid,av,i) )
950	CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, &
951	6, 0 )
952	ENDIF
953
954	END SELECT
955
956	!
957	!-- Select the respective dimension ids
958	IF ( grid_x == 'x' ) THEN
959	id_x = id_dim_x_mask(mid,av)
960	ELSEIF ( grid_x == 'xu' ) THEN
961	id_x = id_dim_xu_mask(mid,av)
962	ENDIF
963
964	IF ( grid_y == 'y' ) THEN
965	id_y = id_dim_y_mask(mid,av)
966	ELSEIF ( grid_y == 'yv' ) THEN
967	id_y = id_dim_yv_mask(mid,av)
968	ENDIF
969
970	IF ( grid_z == 'zu' ) THEN
971	id_z = id_dim_zu_mask(mid,av)
972	ELSEIF ( grid_z == 'zw' ) THEN
973	id_z = id_dim_zw_mask(mid,av)
974	ELSEIF ( grid_z == "zs" ) THEN
975	id_z = id_dim_zs_mask(mid,av)
976	ENDIF
977
978	!
979	!-- Define the grid
980	CALL netcdf_create_var( id_set_mask(mid,av), (/ id_x, id_y, id_z, &
981	id_dim_time_mask(mid,av) /), &
982	domask(mid,av,i), nc_precision(11), &
983	id_var_domask(mid,av,i), &
984	TRIM( domask_unit(mid,av,i) ), &
985	domask(mid,av,i), 494, 495, 496, .TRUE. )
986
987	var_list = TRIM( var_list ) // TRIM( domask(mid,av,i) ) // ';'
988
989	i = i + 1
990
991	ENDDO
992
993	!
994	!-- No arrays to output
995	IF ( i == 1 ) RETURN
996
997	!
998	!-- Write the list of variables as global attribute (this is used by
999	!-- restart runs and by combine_plot_fields)
1000	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
1001	'VAR_LIST', var_list )
1002	CALL netcdf_handle_error( 'netcdf_define_header', 497 )
1003
1004	!
1005	!-- Leave netCDF define mode
1006	nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
1007	CALL netcdf_handle_error( 'netcdf_define_header', 498 )
1008
1009	!
1010	!-- Write data for x (shifted by +dx/2) and xu axis
1011	ALLOCATE( netcdf_data(mask_size(mid,1)) )
1012
1013	netcdf_data = ( mask_i_global(mid,:mask_size(mid,1)) + 0.5_wp ) * dx
1014
1015	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_x_mask(mid,av), &
1016	netcdf_data, start = (/ 1 /), &
1017	count = (/ mask_size(mid,1) /) )
1018	CALL netcdf_handle_error( 'netcdf_define_header', 499 )
1019
1020	netcdf_data = mask_i_global(mid,:mask_size(mid,1)) * dx
1021
1022	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_xu_mask(mid,av),&
1023	netcdf_data, start = (/ 1 /), &
1024	count = (/ mask_size(mid,1) /) )
1025	CALL netcdf_handle_error( 'netcdf_define_header', 500 )
1026
1027	DEALLOCATE( netcdf_data )
1028
1029	!
1030	!-- Write data for y (shifted by +dy/2) and yv axis
1031	ALLOCATE( netcdf_data(mask_size(mid,2)) )
1032
1033	netcdf_data = ( mask_j_global(mid,:mask_size(mid,2)) + 0.5_wp ) * dy
1034
1035	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_y_mask(mid,av), &
1036	netcdf_data, start = (/ 1 /), &
1037	count = (/ mask_size(mid,2) /))
1038	CALL netcdf_handle_error( 'netcdf_define_header', 501 )
1039
1040	netcdf_data = mask_j_global(mid,:mask_size(mid,2)) * dy
1041
1042	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_yv_mask(mid,av), &
1043	netcdf_data, start = (/ 1 /), &
1044	count = (/ mask_size(mid,2) /))
1045	CALL netcdf_handle_error( 'netcdf_define_header', 502 )
1046
1047	DEALLOCATE( netcdf_data )
1048
1049	!
1050	!-- Write zu and zw data (vertical axes)
1051	ALLOCATE( netcdf_data(mask_size(mid,3)) )
1052
1053	netcdf_data = zu( mask_k_global(mid,:mask_size(mid,3)) )
1054
1055	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
1056	netcdf_data, start = (/ 1 /), &
1057	count = (/ mask_size(mid,3) /) )
1058	CALL netcdf_handle_error( 'netcdf_define_header', 503 )
1059
1060	netcdf_data = zw( mask_k_global(mid,:mask_size(mid,3)) )
1061
1062	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
1063	netcdf_data, start = (/ 1 /), &
1064	count = (/ mask_size(mid,3) /) )
1065	CALL netcdf_handle_error( 'netcdf_define_header', 504 )
1066
1067	DEALLOCATE( netcdf_data )
1068
1069	!
1070	!-- In case of non-flat topography write height information
1071	IF ( TRIM( topography ) /= 'flat' .AND. &
1072	netcdf_data_format > 4 ) THEN
1073
1074	ALLOCATE( netcdf_data_2d(mask_size_l(mid,1),mask_size_l(mid,2)) )
1075	netcdf_data_2d = zu_s_inner( mask_i(mid,:mask_size_l(mid,1)), &
1076	mask_j(mid,:mask_size_l(mid,2)) )
1077
1078	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
1079	id_var_zusi_mask(mid,av), &
1080	netcdf_data_2d, &
1081	start = (/ 1, 1 /), &
1082	count = (/ mask_size_l(mid,1), &
1083	mask_size_l(mid,2) /) )
1084	CALL netcdf_handle_error( 'netcdf_define_header', 505 )
1085
1086	netcdf_data_2d = zw_w_inner( mask_i(mid,:mask_size_l(mid,1)), &
1087	mask_j(mid,:mask_size_l(mid,2)) )
1088
1089	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
1090	id_var_zwwi_mask(mid,av), &
1091	netcdf_data_2d, &
1092	start = (/ 1, 1 /), &
1093	count = (/ mask_size_l(mid,1), &
1094	mask_size_l(mid,2) /) )
1095	CALL netcdf_handle_error( 'netcdf_define_header', 506 )
1096
1097	DEALLOCATE( netcdf_data_2d )
1098
1099	ENDIF
1100
1101	IF ( land_surface ) THEN
1102	!
1103	!-- Write zs data (vertical axes for soil model), use negative values
1104	!-- to indicate soil depth
1105	ALLOCATE( netcdf_data(mask_size(mid,3)) )
1106
1107	netcdf_data = zs( mask_k_global(mid,:mask_size(mid,3)) )
1108
1109	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zs_mask(mid,av), &
1110	netcdf_data, start = (/ 1 /), &
1111	count = (/ mask_size(mid,3) /) )
1112	CALL netcdf_handle_error( 'netcdf_define_header', 538 )
1113
1114	DEALLOCATE( netcdf_data )
1115
1116	ENDIF
1117
1118	!
1119	!-- restore original parameter file_id (=formal parameter av) into av
1120	av = file_id
1121
1122
1123	CASE ( 'ma_ext' )
1124
1125	!
1126	!-- decompose actual parameter file_id (=formal parameter av) into
1127	!-- mid and av
1128	file_id = av
1129	IF ( file_id <= 200+max_masks ) THEN
1130	mid = file_id - 200
1131	av = 0
1132	ELSE
1133	mid = file_id - (200+max_masks)
1134	av = 1
1135	ENDIF
1136
1137	!
1138	!-- Get the list of variables and compare with the actual run.
1139	!-- First var_list_old has to be reset, since GET_ATT does not assign
1140	!-- trailing blanks.
1141	var_list_old = ' '
1142	nc_stat = NF90_GET_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'VAR_LIST',&
1143	var_list_old )
1144	CALL netcdf_handle_error( 'netcdf_define_header', 507 )
1145
1146	var_list = ';'
1147	i = 1
1148	DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
1149	var_list = TRIM(var_list) // TRIM( domask(mid,av,i) ) // ';'
1150	i = i + 1
1151	ENDDO
1152
1153	IF ( av == 0 ) THEN
1154	var = '(mask)'
1155	ELSE
1156	var = '(mask_av)'
1157	ENDIF
1158
1159	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
1160	WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ), &
1161	' data for mask', mid, ' from previous run found,', &
1162	'&but this file cannot be extended due to variable ', &
1163	'mismatch.&New file is created instead.'
1164	CALL message( 'define_netcdf_header', 'PA0335', 0, 1, 0, 6, 0 )
1165	extend = .FALSE.
1166	RETURN
1167	ENDIF
1168
1169	!
1170	!-- Get and compare the number of vertical gridpoints
1171	nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'zu_3d', &
1172	id_var_zu_mask(mid,av) )
1173	CALL netcdf_handle_error( 'netcdf_define_header', 508 )
1174
1175	nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av), &
1176	id_var_zu_mask(mid,av), &
1177	dimids = id_dim_zu_mask_old )
1178	CALL netcdf_handle_error( 'netcdf_define_header', 509 )
1179	id_dim_zu_mask(mid,av) = id_dim_zu_mask_old(1)
1180
1181	nc_stat = NF90_INQUIRE_DIMENSION( id_set_mask(mid,av), &
1182	id_dim_zu_mask(mid,av), &
1183	len = nz_old )
1184	CALL netcdf_handle_error( 'netcdf_define_header', 510 )
1185
1186	IF ( mask_size(mid,3) /= nz_old ) THEN
1187	WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ), &
1188	' data for mask', mid, ' from previous run found,', &
1189	'&but this file cannot be extended due to mismatch in ', &
1190	' number of&vertical grid points.', &
1191	'&New file is created instead.'
1192	CALL message( 'define_netcdf_header', 'PA0336', 0, 1, 0, 6, 0 )
1193	extend = .FALSE.
1194	RETURN
1195	ENDIF
1196
1197	!
1198	!-- Get the id of the time coordinate (unlimited coordinate) and its
1199	!-- last index on the file. The next time level is plmask..count+1.
1200	!-- The current time must be larger than the last output time
1201	!-- on the file.
1202	nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'time', &
1203	id_var_time_mask(mid,av) )
1204	CALL netcdf_handle_error( 'netcdf_define_header', 511 )
1205
1206	nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av), &
1207	id_var_time_mask(mid,av), &
1208	dimids = id_dim_time_old )
1209	CALL netcdf_handle_error( 'netcdf_define_header', 512 )
1210	id_dim_time_mask(mid,av) = id_dim_time_old(1)
1211
1212	nc_stat = NF90_INQUIRE_DIMENSION( id_set_mask(mid,av), &
1213	id_dim_time_mask(mid,av), &
1214	len = domask_time_count(mid,av) )
1215	CALL netcdf_handle_error( 'netcdf_define_header', 513 )
1216
1217	nc_stat = NF90_GET_VAR( id_set_mask(mid,av), &
1218	id_var_time_mask(mid,av), &
1219	last_time_coordinate, &
1220	start = (/ domask_time_count(mid,av) /), &
1221	count = (/ 1 /) )
1222	CALL netcdf_handle_error( 'netcdf_define_header', 514 )
1223
1224	IF ( last_time_coordinate(1) >= simulated_time ) 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 because the current ', &
1228	'output time&is less or equal than the last output time ', &
1229	'on this file.&New file is created instead.'
1230	CALL message( 'define_netcdf_header', 'PA0337', 0, 1, 0, 6, 0 )
1231	domask_time_count(mid,av) = 0
1232	extend = .FALSE.
1233	RETURN
1234	ENDIF
1235
1236	!
1237	!-- Dataset seems to be extendable.
1238	!-- Now get the variable ids.
1239	i = 1
1240	DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
1241	nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), &
1242	TRIM( domask(mid,av,i) ), &
1243	id_var_domask(mid,av,i) )
1244	CALL netcdf_handle_error( 'netcdf_define_header', 515 )
1245	i = i + 1
1246	ENDDO
1247
1248	!
1249	!-- Update the title attribute on file
1250	!-- In order to avoid 'data mode' errors if updated attributes are larger
1251	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
1252	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
1253	!-- performance loss due to data copying; an alternative strategy would be
1254	!-- to ensure equal attribute size in a job chain. Maybe revise later.
1255	IF ( av == 0 ) THEN
1256	time_average_text = ' '
1257	ELSE
1258	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
1259	averaging_interval
1260	ENDIF
1261	nc_stat = NF90_REDEF( id_set_mask(mid,av) )
1262	CALL netcdf_handle_error( 'netcdf_define_header', 516 )
1263	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'title', &
1264	TRIM( run_description_header ) // &
1265	TRIM( time_average_text ) )
1266	CALL netcdf_handle_error( 'netcdf_define_header', 517 )
1267	nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
1268	CALL netcdf_handle_error( 'netcdf_define_header', 518 )
1269	WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ), &
1270	' data for mask', mid, ' from previous run found.', &
1271	'&This file will be extended.'
1272	CALL message( 'define_netcdf_header', 'PA0338', 0, 0, 0, 6, 0 )
1273	!
1274	!-- restore original parameter file_id (=formal parameter av) into av
1275	av = file_id
1276
1277
1278	CASE ( '3d_new' )
1279
1280	!
1281	!-- Define some global attributes of the dataset
1282	nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'Conventions', &
1283	'COARDS' )
1284	CALL netcdf_handle_error( 'netcdf_define_header', 62 )
1285
1286	IF ( av == 0 ) THEN
1287	time_average_text = ' '
1288	ELSE
1289	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
1290	averaging_interval
1291	ENDIF
1292	nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'title', &
1293	TRIM( run_description_header ) // &
1294	TRIM( time_average_text ) )
1295	CALL netcdf_handle_error( 'netcdf_define_header', 63 )
1296	IF ( av == 1 ) THEN
1297	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval
1298	nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'time_avg', &
1299	TRIM( time_average_text ) )
1300	CALL netcdf_handle_error( 'netcdf_define_header', 63 )
1301	ENDIF
1302
1303	!
1304	!-- Define time coordinate for volume data.
1305	!-- For parallel output the time dimensions has to be limited, otherwise
1306	!-- the performance drops significantly.
1307	IF ( netcdf_data_format < 5 ) THEN
1308	CALL netcdf_create_dim( id_set_3d(av), 'time', NF90_UNLIMITED, &
1309	id_dim_time_3d(av), 64 )
1310	ELSE
1311	CALL netcdf_create_dim( id_set_3d(av), 'time', ntdim_3d(av), &
1312	id_dim_time_3d(av), 523 )
1313	ENDIF
1314
1315	CALL netcdf_create_var( id_set_3d(av), (/ id_dim_time_3d(av) /), &
1316	'time', NF90_DOUBLE, id_var_time_3d(av), &
1317	'seconds', '', 65, 66, 00 )
1318	!
1319	!-- Define spatial dimensions and coordinates:
1320	!-- Define vertical coordinate grid (zu grid)
1321	CALL netcdf_create_dim( id_set_3d(av), 'zu_3d', nz_do3d-nzb+1, &
1322	id_dim_zu_3d(av), 67 )
1323	CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zu_3d(av) /), &
1324	'zu_3d', NF90_DOUBLE, id_var_zu_3d(av), &
1325	'meters', '', 68, 69, 00 )
1326	!
1327	!-- Define vertical coordinate grid (zw grid)
1328	CALL netcdf_create_dim( id_set_3d(av), 'zw_3d', nz_do3d-nzb+1, &
1329	id_dim_zw_3d(av), 70 )
1330	CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zw_3d(av) /), &
1331	'zw_3d', NF90_DOUBLE, id_var_zw_3d(av), &
1332	'meters', '', 71, 72, 00 )
1333	!
1334	!-- Define x-axis (for scalar position)
1335	CALL netcdf_create_dim( id_set_3d(av), 'x', nx+1, id_dim_x_3d(av), &
1336	73 )
1337	CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av) /), 'x', &
1338	NF90_DOUBLE, id_var_x_3d(av), 'meters', '', &
1339	74, 75, 00 )
1340	!
1341	!-- Define x-axis (for u position)
1342	CALL netcdf_create_dim( id_set_3d(av), 'xu', nx+1, id_dim_xu_3d(av), &
1343	358 )
1344	CALL netcdf_create_var( id_set_3d(av), (/ id_dim_xu_3d(av) /), 'xu', &
1345	NF90_DOUBLE, id_var_xu_3d(av), 'meters', '', &
1346	359, 360, 000 )
1347	!
1348	!-- Define y-axis (for scalar position)
1349	CALL netcdf_create_dim( id_set_3d(av), 'y', ny+1, id_dim_y_3d(av), &
1350	76 )
1351	CALL netcdf_create_var( id_set_3d(av), (/ id_dim_y_3d(av) /), 'y', &
1352	NF90_DOUBLE, id_var_y_3d(av), 'meters', '', &
1353	77, 78, 00 )
1354	!
1355	!-- Define y-axis (for v position)
1356	CALL netcdf_create_dim( id_set_3d(av), 'yv', ny+1, id_dim_yv_3d(av), &
1357	361 )
1358	CALL netcdf_create_var( id_set_3d(av), (/ id_dim_yv_3d(av) /), 'yv', &
1359	NF90_DOUBLE, id_var_yv_3d(av), 'meters', '', &
1360	362, 363, 000 )
1361	!
1362	!-- In case of non-flat topography define 2d-arrays containing the height
1363	!-- information. Only output 2d topography information in case of parallel
1364	!-- output.
1365	IF ( TRIM( topography ) /= 'flat' .AND. &
1366	netcdf_data_format > 4 ) THEN
1367	!
1368	!-- Define zusi = zu(nzb_s_inner)
1369	CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av), &
1370	id_dim_y_3d(av) /), 'zusi', NF90_DOUBLE, &
1371	id_var_zusi_3d(av), 'meters', &
1372	'zu(nzb_s_inner)', 413, 414, 415 )
1373	!
1374	!-- Define zwwi = zw(nzb_w_inner)
1375	CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av), &
1376	id_dim_y_3d(av) /), 'zwwi', NF90_DOUBLE, &
1377	id_var_zwwi_3d(av), 'meters', &
1378	'zw(nzb_w_inner)', 416, 417, 418 )
1379
1380	ENDIF
1381
1382	IF ( land_surface ) THEN
1383	!
1384	!-- Define vertical coordinate grid (zs grid)
1385	CALL netcdf_create_dim( id_set_3d(av), 'zs_3d', &
1386	nzt_soil-nzb_soil+1, id_dim_zs_3d(av), 70 )
1387	CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zs_3d(av) /), &
1388	'zs_3d', NF90_DOUBLE, id_var_zs_3d(av), &
1389	'meters', '', 71, 72, 00 )
1390
1391	ENDIF
1392
1393	!
1394	!-- Define the variables
1395	var_list = ';'
1396	i = 1
1397
1398	DO WHILE ( do3d(av,i)(1:1) /= ' ' )
1399	!
1400	!-- Temporary solution to account for data output within the new urban
1401	!-- surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
1402	trimvar = TRIM( do3d(av,i) )
1403	IF ( urban_surface .AND. trimvar(1:4) == 'usm_' ) THEN
1404	trimvar = 'usm_output'
1405	ENDIF
1406	!
1407	!-- Check for the grid
1408	found = .FALSE.
1409	SELECT CASE ( trimvar )
1410	!
1411	!-- Most variables are defined on the scalar grid
1412	CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', &
1413	'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', &
1414	'rho_ocean', 's', 'sa', 'vpt' )
1415
1416	grid_x = 'x'
1417	grid_y = 'y'
1418	grid_z = 'zu'
1419	!
1420	!-- u grid
1421	CASE ( 'u' )
1422
1423	grid_x = 'xu'
1424	grid_y = 'y'
1425	grid_z = 'zu'
1426	!
1427	!-- v grid
1428	CASE ( 'v' )
1429
1430	grid_x = 'x'
1431	grid_y = 'yv'
1432	grid_z = 'zu'
1433	!
1434	!-- w grid
1435	CASE ( 'w' )
1436
1437	grid_x = 'x'
1438	grid_y = 'y'
1439	grid_z = 'zw'
1440
1441	!
1442	!-- Block of urban surface model outputs
1443	CASE ( 'usm_output' )
1444	CALL usm_define_netcdf_grid( do3d(av,i), found, &
1445	grid_x, grid_y, grid_z )
1446
1447	CASE DEFAULT
1448
1449	CALL tcm_define_netcdf_grid( do3d(av,i), found, &
1450	grid_x, grid_y, grid_z )
1451
1452	!
1453	!-- Check for land surface quantities
1454	IF ( land_surface ) THEN
1455	CALL lsm_define_netcdf_grid( do3d(av,i), found, grid_x, &
1456	grid_y, grid_z )
1457	ENDIF
1458
1459	!
1460	!-- Check for plant canopy quantities
1461	IF ( .NOT. found .AND. plant_canopy ) THEN
1462	CALL pcm_define_netcdf_grid( do3d(av,i), found, grid_x, &
1463	grid_y, grid_z )
1464	ENDIF
1465
1466	!
1467	!-- Check for radiation quantities
1468	IF ( .NOT. found .AND. radiation ) THEN
1469	CALL radiation_define_netcdf_grid( do3d(av,i), found, &
1470	grid_x, grid_y, &
1471	grid_z )
1472	ENDIF
1473
1474	!
1475	!-- Check for chemistry quantities
1476	#if defined( __chem )
1477	IF ( .NOT. found .AND. air_chemistry ) THEN
1478	CALL chem_define_netcdf_grid( do3d(av,i), found, &
1479	grid_x, grid_y, grid_z )
1480	ENDIF
1481	#endif
1482
1483	!-- Check for user-defined quantities
1484	IF ( .NOT. found ) THEN
1485	CALL user_define_netcdf_grid( do3d(av,i), found, grid_x, &
1486	grid_y, grid_z )
1487	ENDIF
1488
1489	IF ( .NOT. found ) THEN
1490	WRITE ( message_string, * ) 'no grid defined for varia', &
1491	'ble ', TRIM( do3d(av,i) )
1492	CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, &
1493	6, 0 )
1494	ENDIF
1495
1496	END SELECT
1497
1498	!
1499	!-- Select the respective dimension ids
1500	IF ( grid_x == 'x' ) THEN
1501	id_x = id_dim_x_3d(av)
1502	ELSEIF ( grid_x == 'xu' ) THEN
1503	id_x = id_dim_xu_3d(av)
1504	ENDIF
1505
1506	IF ( grid_y == 'y' ) THEN
1507	id_y = id_dim_y_3d(av)
1508	ELSEIF ( grid_y == 'yv' ) THEN
1509	id_y = id_dim_yv_3d(av)
1510	ENDIF
1511
1512	IF ( grid_z == 'zu' ) THEN
1513	id_z = id_dim_zu_3d(av)
1514	ELSEIF ( grid_z == 'zw' ) THEN
1515	id_z = id_dim_zw_3d(av)
1516	ELSEIF ( grid_z == 'zs' ) THEN
1517	id_z = id_dim_zs_3d(av)
1518	ENDIF
1519
1520	!
1521	!-- Define the grid
1522	CALL netcdf_create_var( id_set_3d(av),(/ id_x, id_y, id_z, &
1523	id_dim_time_3d(av) /), do3d(av,i), &
1524	nc_precision(4), id_var_do3d(av,i), &
1525	TRIM( do3d_unit(av,i) ), do3d(av,i), 79, &
1526	80, 357, .TRUE. )
1527	#if defined( __netcdf4_parallel )
1528	IF ( netcdf_data_format > 4 ) THEN
1529	!
1530	!-- Set no fill for every variable to increase performance.
1531	nc_stat = NF90_DEF_VAR_FILL( id_set_3d(av), &
1532	id_var_do3d(av,i), &
1533	1, 0 )
1534	CALL netcdf_handle_error( 'netcdf_define_header', 532 )
1535	!
1536	!-- Set collective io operations for parallel io
1537	nc_stat = NF90_VAR_PAR_ACCESS( id_set_3d(av), &
1538	id_var_do3d(av,i), &
1539	NF90_COLLECTIVE )
1540	CALL netcdf_handle_error( 'netcdf_define_header', 445 )
1541	ENDIF
1542	#endif
1543	var_list = TRIM( var_list ) // TRIM( do3d(av,i) ) // ';'
1544
1545	i = i + 1
1546
1547	ENDDO
1548
1549	!
1550	!-- No arrays to output
1551	IF ( i == 1 ) RETURN
1552
1553	!
1554	!-- Write the list of variables as global attribute (this is used by
1555	!-- restart runs and by combine_plot_fields)
1556	nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
1557	var_list )
1558	CALL netcdf_handle_error( 'netcdf_define_header', 81 )
1559
1560	!
1561	!-- Set general no fill, otherwise the performance drops significantly for
1562	!-- parallel output.
1563	nc_stat = NF90_SET_FILL( id_set_3d(av), NF90_NOFILL, oldmode )
1564	CALL netcdf_handle_error( 'netcdf_define_header', 528 )
1565
1566	!
1567	!-- Leave netCDF define mode
1568	nc_stat = NF90_ENDDEF( id_set_3d(av) )
1569	CALL netcdf_handle_error( 'netcdf_define_header', 82 )
1570
1571	!
1572	!-- These data are only written by PE0 for parallel output to increase
1573	!-- the performance.
1574	IF ( myid == 0 .OR. netcdf_data_format < 5 ) THEN
1575	!
1576	!-- Write data for x (shifted by +dx/2) and xu axis
1577	ALLOCATE( netcdf_data(0:nx) )
1578
1579	DO i = 0, nx
1580	netcdf_data(i) = ( i + 0.5 ) * dx
1581	ENDDO
1582
1583	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_x_3d(av), &
1584	netcdf_data, start = (/ 1 /), &
1585	count = (/ nx+1 /) )
1586	CALL netcdf_handle_error( 'netcdf_define_header', 83 )
1587
1588	DO i = 0, nx
1589	netcdf_data(i) = i * dx
1590	ENDDO
1591
1592	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_xu_3d(av), &
1593	netcdf_data, start = (/ 1 /), &
1594	count = (/ nx+1 /) )
1595	CALL netcdf_handle_error( 'netcdf_define_header', 385 )
1596
1597	DEALLOCATE( netcdf_data )
1598
1599	!
1600	!-- Write data for y (shifted by +dy/2) and yv axis
1601	ALLOCATE( netcdf_data(0:ny) )
1602
1603	DO i = 0, ny
1604	netcdf_data(i) = ( i + 0.5_wp ) * dy
1605	ENDDO
1606
1607	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_y_3d(av), &
1608	netcdf_data, start = (/ 1 /), &
1609	count = (/ ny+1 /) )
1610	CALL netcdf_handle_error( 'netcdf_define_header', 84 )
1611
1612	DO i = 0, ny
1613	netcdf_data(i) = i * dy
1614	ENDDO
1615
1616	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_yv_3d(av), &
1617	netcdf_data, start = (/ 1 /), &
1618	count = (/ ny+1 /))
1619	CALL netcdf_handle_error( 'netcdf_define_header', 387 )
1620
1621	DEALLOCATE( netcdf_data )
1622
1623	!
1624	!-- Write zu and zw data (vertical axes)
1625	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zu_3d(av), &
1626	zu(nzb:nz_do3d), start = (/ 1 /), &
1627	count = (/ nz_do3d-nzb+1 /) )
1628	CALL netcdf_handle_error( 'netcdf_define_header', 85 )
1629
1630
1631	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zw_3d(av), &
1632	zw(nzb:nz_do3d), start = (/ 1 /), &
1633	count = (/ nz_do3d-nzb+1 /) )
1634	CALL netcdf_handle_error( 'netcdf_define_header', 86 )
1635
1636	IF ( land_surface ) THEN
1637	!
1638	!-- Write zs grid
1639	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zs_3d(av), &
1640	- zs(nzb_soil:nzt_soil), start = (/ 1 /), &
1641	count = (/ nzt_soil-nzb_soil+1 /) )
1642	CALL netcdf_handle_error( 'netcdf_define_header', 86 )
1643	ENDIF
1644
1645	ENDIF
1646	!
1647	!-- In case of non-flat topography write height information. Only for
1648	!-- parallel netcdf output.
1649	IF ( TRIM( topography ) /= 'flat' .AND. &
1650	netcdf_data_format > 4 ) THEN
1651
1652	! IF ( nxr == nx .AND. nyn /= ny ) THEN
1653	! nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av), &
1654	! zu_s_inner(nxl:nxr+1,nys:nyn), &
1655	! start = (/ nxl+1, nys+1 /), &
1656	! count = (/ nxr-nxl+2, nyn-nys+1 /) )
1657	! ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN
1658	! nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av), &
1659	! zu_s_inner(nxl:nxr,nys:nyn+1), &
1660	! start = (/ nxl+1, nys+1 /), &
1661	! count = (/ nxr-nxl+1, nyn-nys+2 /) )
1662	! ELSEIF ( nxr == nx .AND. nyn == ny ) THEN
1663	! nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av), &
1664	! zu_s_inner(nxl:nxr+1,nys:nyn+1), &
1665	! start = (/ nxl+1, nys+1 /), &
1666	! count = (/ nxr-nxl+2, nyn-nys+2 /) )
1667	! ELSE
1668	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av), &
1669	zu_s_inner(nxl:nxr,nys:nyn), &
1670	start = (/ nxl+1, nys+1 /), &
1671	count = (/ nxr-nxl+1, nyn-nys+1 /) )
1672	! ENDIF
1673	CALL netcdf_handle_error( 'netcdf_define_header', 419 )
1674
1675	! IF ( nxr == nx .AND. nyn /= ny ) THEN
1676	! nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av), &
1677	! zw_w_inner(nxl:nxr+1,nys:nyn), &
1678	! start = (/ nxl+1, nys+1 /), &
1679	! count = (/ nxr-nxl+2, nyn-nys+1 /) )
1680	! ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN
1681	! nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av), &
1682	! zw_w_inner(nxl:nxr,nys:nyn+1), &
1683	! start = (/ nxl+1, nys+1 /), &
1684	! count = (/ nxr-nxl+1, nyn-nys+2 /) )
1685	! ELSEIF ( nxr == nx .AND. nyn == ny ) THEN
1686	! nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av), &
1687	! zw_w_inner(nxl:nxr+1,nys:nyn+1), &
1688	! start = (/ nxl+1, nys+1 /), &
1689	! count = (/ nxr-nxl+2, nyn-nys+2 /) )
1690	! ELSE
1691	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av), &
1692	zw_w_inner(nxl:nxr,nys:nyn), &
1693	start = (/ nxl+1, nys+1 /), &
1694	count = (/ nxr-nxl+1, nyn-nys+1 /) )
1695	! ENDIF
1696	CALL netcdf_handle_error( 'netcdf_define_header', 420 )
1697
1698	ENDIF
1699
1700	CASE ( '3d_ext' )
1701
1702	!
1703	!-- Get the list of variables and compare with the actual run.
1704	!-- First var_list_old has to be reset, since GET_ATT does not assign
1705	!-- trailing blanks.
1706	var_list_old = ' '
1707	nc_stat = NF90_GET_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
1708	var_list_old )
1709	CALL netcdf_handle_error( 'netcdf_define_header', 87 )
1710
1711	var_list = ';'
1712	i = 1
1713	DO WHILE ( do3d(av,i)(1:1) /= ' ' )
1714	var_list = TRIM(var_list) // TRIM( do3d(av,i) ) // ';'
1715	i = i + 1
1716	ENDDO
1717
1718	IF ( av == 0 ) THEN
1719	var = '(3d)'
1720	ELSE
1721	var = '(3d_av)'
1722	ENDIF
1723
1724	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
1725	message_string = 'netCDF file for volume data ' // &
1726	TRIM( var ) // ' from previous run found,' // &
1727	'&but this file cannot be extended due to' // &
1728	' variable mismatch.' // &
1729	'&New file is created instead.'
1730	CALL message( 'define_netcdf_header', 'PA0245', 0, 1, 0, 6, 0 )
1731	extend = .FALSE.
1732	RETURN
1733	ENDIF
1734
1735	!
1736	!-- Get and compare the number of vertical gridpoints
1737	nc_stat = NF90_INQ_VARID( id_set_3d(av), 'zu_3d', id_var_zu_3d(av) )
1738	CALL netcdf_handle_error( 'netcdf_define_header', 88 )
1739
1740	nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_zu_3d(av), &
1741	dimids = id_dim_zu_3d_old )
1742	CALL netcdf_handle_error( 'netcdf_define_header', 89 )
1743	id_dim_zu_3d(av) = id_dim_zu_3d_old(1)
1744
1745	nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_zu_3d(av), &
1746	len = nz_old )
1747	CALL netcdf_handle_error( 'netcdf_define_header', 90 )
1748
1749	IF ( nz_do3d-nzb+1 /= nz_old ) THEN
1750	message_string = 'netCDF file for volume data ' // &
1751	TRIM( var ) // ' from previous run found,' // &
1752	'&but this file cannot be extended due to' // &
1753	' mismatch in number of' // &
1754	'&vertical grid points (nz_do3d).' // &
1755	'&New file is created instead.'
1756	CALL message( 'define_netcdf_header', 'PA0246', 0, 1, 0, 6, 0 )
1757	extend = .FALSE.
1758	RETURN
1759	ENDIF
1760
1761	!
1762	!-- Get the id of the time coordinate (unlimited coordinate) and its
1763	!-- last index on the file. The next time level is pl3d..count+1.
1764	!-- The current time must be larger than the last output time
1765	!-- on the file.
1766	nc_stat = NF90_INQ_VARID( id_set_3d(av), 'time', id_var_time_3d(av) )
1767	CALL netcdf_handle_error( 'netcdf_define_header', 91 )
1768
1769	nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_time_3d(av), &
1770	dimids = id_dim_time_old )
1771	CALL netcdf_handle_error( 'netcdf_define_header', 92 )
1772
1773	id_dim_time_3d(av) = id_dim_time_old(1)
1774
1775	nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_time_3d(av), &
1776	len = ntime_count )
1777	CALL netcdf_handle_error( 'netcdf_define_header', 93 )
1778
1779	!
1780	!-- For non-parallel output use the last output time level of the netcdf
1781	!-- file because the time dimension is unlimited. In case of parallel
1782	!-- output the variable ntime_count could get the value of 9*10E36 because
1783	!-- the time dimension is limited.
1784	IF ( netcdf_data_format < 5 ) do3d_time_count(av) = ntime_count
1785
1786	nc_stat = NF90_GET_VAR( id_set_3d(av), id_var_time_3d(av), &
1787	last_time_coordinate, &
1788	start = (/ do3d_time_count(av) /), &
1789	count = (/ 1 /) )
1790	CALL netcdf_handle_error( 'netcdf_define_header', 94 )
1791
1792	IF ( last_time_coordinate(1) >= simulated_time ) THEN
1793	message_string = 'netCDF file for volume data ' // &
1794	TRIM( var ) // ' from previous run found,' // &
1795	'&but this file cannot be extended becaus' // &
1796	'e the current output time' // &
1797	'&is less or equal than the last output t' // &
1798	'ime on this file.' // &
1799	'&New file is created instead.'
1800	CALL message( 'define_netcdf_header', 'PA0247', 0, 1, 0, 6, 0 )
1801	do3d_time_count(av) = 0
1802	extend = .FALSE.
1803	RETURN
1804	ENDIF
1805
1806	IF ( netcdf_data_format > 4 ) THEN
1807	!
1808	!-- Set needed time levels (ntdim) to
1809	!-- saved time levels + to be saved time levels.
1810	IF ( av == 0 ) THEN
1811	ntdim_3d(0) = do3d_time_count(0) + CEILING( &
1812	( end_time - MAX( skip_time_do3d, &
1813	simulated_time_at_begin ) &
1814	) / dt_do3d )
1815	ELSE
1816	ntdim_3d(1) = do3d_time_count(1) + CEILING( &
1817	( end_time - MAX( skip_time_data_output_av, &
1818	simulated_time_at_begin ) &
1819	) / dt_data_output_av )
1820	ENDIF
1821
1822	!
1823	!-- Check if the needed number of output time levels is increased
1824	!-- compared to the number of time levels in the existing file.
1825	IF ( ntdim_3d(av) > ntime_count ) THEN
1826	message_string = 'netCDF file for volume data ' // &
1827	TRIM( var ) // ' from previous run found,' // &
1828	'&but this file cannot be extended becaus' // &
1829	'e the number of output time levels&has b' // &
1830	'een increased compared to the previous s' // &
1831	'imulation.' // &
1832	'&New file is created instead.'
1833	CALL message( 'define_netcdf_header', 'PA0388', 0, 1, 0, 6, 0 )
1834	do3d_time_count(av) = 0
1835	extend = .FALSE.
1836	!
1837	!-- Recalculate the needed time levels for the new file.
1838	IF ( av == 0 ) THEN
1839	ntdim_3d(0) = CEILING( &
1840	( end_time - MAX( skip_time_do3d, &
1841	simulated_time_at_begin ) &
1842	) / dt_do3d )
1843	ELSE
1844	ntdim_3d(1) = CEILING( &
1845	( end_time - MAX( skip_time_data_output_av, &
1846	simulated_time_at_begin ) &
1847	) / dt_data_output_av )
1848	ENDIF
1849	RETURN
1850	ENDIF
1851	ENDIF
1852
1853	!
1854	!-- Dataset seems to be extendable.
1855	!-- Now get the variable ids.
1856	i = 1
1857	DO WHILE ( do3d(av,i)(1:1) /= ' ' )
1858	nc_stat = NF90_INQ_VARID( id_set_3d(av), TRIM( do3d(av,i) ), &
1859	id_var_do3d(av,i) )
1860	CALL netcdf_handle_error( 'netcdf_define_header', 95 )
1861	#if defined( __netcdf4_parallel )
1862	!
1863	!-- Set collective io operations for parallel io
1864	IF ( netcdf_data_format > 4 ) THEN
1865	nc_stat = NF90_VAR_PAR_ACCESS( id_set_3d(av), &
1866	id_var_do3d(av,i), &
1867	NF90_COLLECTIVE )
1868	CALL netcdf_handle_error( 'netcdf_define_header', 453 )
1869	ENDIF
1870	#endif
1871	i = i + 1
1872	ENDDO
1873
1874	!
1875	!-- Update the title attribute on file
1876	!-- In order to avoid 'data mode' errors if updated attributes are larger
1877	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
1878	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
1879	!-- performance loss due to data copying; an alternative strategy would be
1880	!-- to ensure equal attribute size. Maybe revise later.
1881	IF ( av == 0 ) THEN
1882	time_average_text = ' '
1883	ELSE
1884	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
1885	averaging_interval
1886	ENDIF
1887	nc_stat = NF90_REDEF( id_set_3d(av) )
1888	CALL netcdf_handle_error( 'netcdf_define_header', 429 )
1889	nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'title', &
1890	TRIM( run_description_header ) // &
1891	TRIM( time_average_text ) )
1892	CALL netcdf_handle_error( 'netcdf_define_header', 96 )
1893	nc_stat = NF90_ENDDEF( id_set_3d(av) )
1894	CALL netcdf_handle_error( 'netcdf_define_header', 430 )
1895	message_string = 'netCDF file for volume data ' // &
1896	TRIM( var ) // ' from previous run found.' // &
1897	'&This file will be extended.'
1898	CALL message( 'define_netcdf_header', 'PA0248', 0, 0, 0, 6, 0 )
1899
1900	CASE ( 'xy_new' )
1901
1902	!
1903	!-- Define some global attributes of the dataset
1904	nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'Conventions', &
1905	'COARDS' )
1906	CALL netcdf_handle_error( 'netcdf_define_header', 97 )
1907
1908	IF ( av == 0 ) THEN
1909	time_average_text = ' '
1910	ELSE
1911	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
1912	averaging_interval
1913	ENDIF
1914	nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'title', &
1915	TRIM( run_description_header ) // &
1916	TRIM( time_average_text ) )
1917	CALL netcdf_handle_error( 'netcdf_define_header', 98 )
1918	IF ( av == 1 ) THEN
1919	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval
1920	nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'time_avg', &
1921	TRIM( time_average_text ) )
1922	CALL netcdf_handle_error( 'netcdf_define_header', 98 )
1923	ENDIF
1924
1925	!
1926	!-- Define time coordinate for xy sections.
1927	!-- For parallel output the time dimensions has to be limited, otherwise
1928	!-- the performance drops significantly.
1929	IF ( netcdf_data_format < 5 ) THEN
1930	CALL netcdf_create_dim( id_set_xy(av), 'time', NF90_UNLIMITED, &
1931	id_dim_time_xy(av), 99 )
1932	ELSE
1933	CALL netcdf_create_dim( id_set_xy(av), 'time', ntdim_2d_xy(av), &
1934	id_dim_time_xy(av), 524 )
1935	ENDIF
1936
1937	CALL netcdf_create_var( id_set_xy(av), (/ id_dim_time_xy(av) /), &
1938	'time', NF90_DOUBLE, id_var_time_xy(av), &
1939	'seconds', '', 100, 101, 000 )
1940	!
1941	!-- Define the spatial dimensions and coordinates for xy-sections.
1942	!-- First, determine the number of horizontal sections to be written.
1943	IF ( section(1,1) == -9999 ) THEN
1944	RETURN
1945	ELSE
1946	ns = 1
1947	DO WHILE ( section(ns,1) /= -9999 .AND. ns <= 100 )
1948	ns = ns + 1
1949	ENDDO
1950	ns = ns - 1
1951	ENDIF
1952
1953	!
1954	!-- Define vertical coordinate grid (zu grid)
1955	CALL netcdf_create_dim( id_set_xy(av), 'zu_xy', ns, &
1956	id_dim_zu_xy(av), 102 )
1957	CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu_xy(av) /), &
1958	'zu_xy', NF90_DOUBLE, id_var_zu_xy(av), &
1959	'meters', '', 103, 104, 000 )
1960	!
1961	!-- Define vertical coordinate grid (zw grid)
1962	CALL netcdf_create_dim( id_set_xy(av), 'zw_xy', ns, &
1963	id_dim_zw_xy(av), 105 )
1964	CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zw_xy(av) /), &
1965	'zw_xy', NF90_DOUBLE, id_var_zw_xy(av), &
1966	'meters', '', 106, 107, 000 )
1967
1968	IF ( land_surface ) THEN
1969
1970	ns_do = 1
1971	DO WHILE ( section(ns_do,1) /= -9999 .AND. ns_do < nzs )
1972	ns_do = ns_do + 1
1973	ENDDO
1974	!
1975	!-- Define vertical coordinate grid (zs grid)
1976	CALL netcdf_create_dim( id_set_xy(av), 'zs_xy', ns_do, &
1977	id_dim_zs_xy(av), 539 )
1978	CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zs_xy(av) /), &
1979	'zs_xy', NF90_DOUBLE, id_var_zs_xy(av), &
1980	'meters', '', 540, 541, 000 )
1981
1982	ENDIF
1983
1984	!
1985	!-- Define a pseudo vertical coordinate grid for the surface variables
1986	!-- u* and t* to store their height level
1987	CALL netcdf_create_dim( id_set_xy(av), 'zu1_xy', 1, &
1988	id_dim_zu1_xy(av), 108 )
1989	CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu1_xy(av) /), &
1990	'zu1_xy', NF90_DOUBLE, id_var_zu1_xy(av), &
1991	'meters', '', 109, 110, 000 )
1992	!
1993	!-- Define a variable to store the layer indices of the horizontal cross
1994	!-- sections, too
1995	CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu_xy(av) /), &
1996	'ind_z_xy', NF90_DOUBLE, &
1997	id_var_ind_z_xy(av), 'gridpoints', '', 111, &
1998	112, 000 )
1999	!
2000	!-- Define x-axis (for scalar position)
2001	CALL netcdf_create_dim( id_set_xy(av), 'x', nx+1, id_dim_x_xy(av), &
2002	113 )
2003	CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av) /), 'x', &
2004	NF90_DOUBLE, id_var_x_xy(av), 'meters', '', &
2005	114, 115, 000 )
2006	!
2007	!-- Define x-axis (for u position)
2008	CALL netcdf_create_dim( id_set_xy(av), 'xu', nx+1, &
2009	id_dim_xu_xy(av), 388 )
2010	CALL netcdf_create_var( id_set_xy(av), (/ id_dim_xu_xy(av) /), 'xu', &
2011	NF90_DOUBLE, id_var_xu_xy(av), 'meters', '', &
2012	389, 390, 000 )
2013	!
2014	!-- Define y-axis (for scalar position)
2015	CALL netcdf_create_dim( id_set_xy(av), 'y', ny+1, id_dim_y_xy(av), &
2016	116 )
2017	CALL netcdf_create_var( id_set_xy(av), (/ id_dim_y_xy(av) /), 'y', &
2018	NF90_DOUBLE, id_var_y_xy(av), 'meters', '', &
2019	117, 118, 000 )
2020	!
2021	!-- Define y-axis (for scalar position)
2022	CALL netcdf_create_dim( id_set_xy(av), 'yv', ny+1, &
2023	id_dim_yv_xy(av), 364 )
2024	CALL netcdf_create_var( id_set_xy(av), (/ id_dim_yv_xy(av) /), 'yv', &
2025	NF90_DOUBLE, id_var_yv_xy(av), 'meters', '', &
2026	365, 366, 000 )
2027	!
2028	!-- In case of non-flat topography define 2d-arrays containing the height
2029	!-- information. Only for parallel netcdf output.
2030	IF ( TRIM( topography ) /= 'flat' .AND. &
2031	netcdf_data_format > 4 ) THEN
2032	!
2033	!-- Define zusi = zu(nzb_s_inner)
2034	CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av), &
2035	id_dim_y_xy(av) /), 'zusi', NF90_DOUBLE, &
2036	id_var_zusi_xy(av), 'meters', &
2037	'zu(nzb_s_inner)', 421, 422, 423 )
2038	!
2039	!-- Define zwwi = zw(nzb_w_inner)
2040	CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av), &
2041	id_dim_y_xy(av) /), 'zwwi', NF90_DOUBLE, &
2042	id_var_zwwi_xy(av), 'meters', &
2043	'zw(nzb_w_inner)', 424, 425, 426 )
2044
2045	ENDIF
2046
2047	!
2048	!-- Define the variables
2049	var_list = ';'
2050	i = 1
2051
2052	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2053
2054	IF ( INDEX( do2d(av,i), 'xy' ) /= 0 ) THEN
2055	!
2056	!-- If there is a star in the variable name (u* or t*), it is a
2057	!-- surface variable. Define it with id_dim_zu1_xy.
2058	IF ( INDEX( do2d(av,i), '*' ) /= 0 ) THEN
2059
2060	CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av), &
2061	id_dim_y_xy(av), id_dim_zu1_xy(av), &
2062	id_dim_time_xy(av) /), do2d(av,i), &
2063	nc_precision(1), id_var_do2d(av,i), &
2064	TRIM( do2d_unit(av,i) ), &
2065	do2d(av,i), 119, 120, 354, .TRUE. )
2066
2067	ELSE
2068
2069	!
2070	!-- Check for the grid
2071	found = .FALSE.
2072	SELECT CASE ( do2d(av,i) )
2073	!
2074	!-- Most variables are defined on the zu grid
2075	CASE ( 'e_xy', 'lpt_xy', 'nc_xy', 'nr_xy', 'p_xy', &
2076	'pc_xy', 'pr_xy', 'prr_xy', 'pt_xy', 'q_xy', &
2077	'qc_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy', &
2078	'ql_vp_xy', 'qr_xy', 'qv_xy', 'rho_ocean_xy', &
2079	's_xy', 'sa_xy', 'vpt_xy' )
2080
2081	grid_x = 'x'
2082	grid_y = 'y'
2083	grid_z = 'zu'
2084	!
2085	!-- u grid
2086	CASE ( 'u_xy' )
2087
2088	grid_x = 'xu'
2089	grid_y = 'y'
2090	grid_z = 'zu'
2091	!
2092	!-- v grid
2093	CASE ( 'v_xy' )
2094
2095	grid_x = 'x'
2096	grid_y = 'yv'
2097	grid_z = 'zu'
2098	!
2099	!-- w grid
2100	CASE ( 'w_xy' )
2101
2102	grid_x = 'x'
2103	grid_y = 'y'
2104	grid_z = 'zw'
2105
2106
2107	CASE DEFAULT
2108	!
2109	!-- Check for land surface quantities
2110	IF ( land_surface ) THEN
2111	CALL lsm_define_netcdf_grid( do2d(av,i), found, &
2112	grid_x, grid_y, grid_z )
2113	ENDIF
2114
2115	IF ( .NOT. found ) THEN
2116	CALL tcm_define_netcdf_grid( do2d(av,i), found, &
2117	grid_x, grid_y, &
2118	grid_z )
2119	ENDIF
2120
2121	!
2122	!-- Check for radiation quantities
2123	IF ( .NOT. found .AND. radiation ) THEN
2124	CALL radiation_define_netcdf_grid( do2d(av,i), &
2125	found, grid_x, grid_y,&
2126	grid_z )
2127	ENDIF
2128
2129	!
2130	!-- Check for chemistry quantities
2131	#if defined( __chem )
2132	IF ( .NOT. found .AND. air_chemistry ) THEN
2133	CALL chem_define_netcdf_grid( do2d(av,i), found, &
2134	grid_x, grid_y, &
2135	grid_z )
2136	ENDIF
2137	#endif
2138
2139	!
2140	!-- Check for UV exposure quantities
2141	IF ( .NOT. found .AND. uv_exposure ) THEN
2142	CALL uvem_define_netcdf_grid( do2d(av,i), found, &
2143	grid_x, grid_y, grid_z )
2144	ENDIF
2145
2146	!
2147	!-- Check for user-defined quantities
2148	IF ( .NOT. found ) THEN
2149	CALL user_define_netcdf_grid( do2d(av,i), found, &
2150	grid_x, grid_y, &
2151	grid_z )
2152	ENDIF
2153
2154	IF ( .NOT. found ) THEN
2155	WRITE ( message_string, * ) 'no grid defined for', &
2156	' variable ', TRIM( do2d(av,i) )
2157	CALL message( 'define_netcdf_header', 'PA0244', &
2158	0, 1, 0, 6, 0 )
2159	ENDIF
2160
2161	END SELECT
2162
2163	!
2164	!-- Select the respective dimension ids
2165	IF ( grid_x == 'x' ) THEN
2166	id_x = id_dim_x_xy(av)
2167	ELSEIF ( grid_x == 'xu' ) THEN
2168	id_x = id_dim_xu_xy(av)
2169	ENDIF
2170
2171	IF ( grid_y == 'y' ) THEN
2172	id_y = id_dim_y_xy(av)
2173	ELSEIF ( grid_y == 'yv' ) THEN
2174	id_y = id_dim_yv_xy(av)
2175	ENDIF
2176
2177	IF ( grid_z == 'zu' ) THEN
2178	id_z = id_dim_zu_xy(av)
2179	ELSEIF ( grid_z == 'zw' ) THEN
2180	id_z = id_dim_zw_xy(av)
2181	ELSEIF ( grid_z == 'zs' ) THEN
2182	id_z = id_dim_zs_xy(av)
2183	ENDIF
2184
2185	!
2186	!-- Define the grid
2187	CALL netcdf_create_var( id_set_xy(av), (/ id_x, id_y, id_z, &
2188	id_dim_time_xy(av) /), do2d(av,i), &
2189	nc_precision(1), id_var_do2d(av,i), &
2190	TRIM( do2d_unit(av,i) ), &
2191	do2d(av,i), 119, 120, 354, .TRUE. )
2192
2193	ENDIF
2194
2195	#if defined( __netcdf4_parallel )
2196	IF ( netcdf_data_format > 4 ) THEN
2197	!
2198	!-- Set no fill for every variable to increase performance.
2199	nc_stat = NF90_DEF_VAR_FILL( id_set_xy(av), &
2200	id_var_do2d(av,i), &
2201	1, 0 )
2202	CALL netcdf_handle_error( 'netcdf_define_header', 533 )
2203	!
2204	!-- Set collective io operations for parallel io
2205	nc_stat = NF90_VAR_PAR_ACCESS( id_set_xy(av), &
2206	id_var_do2d(av,i), &
2207	NF90_COLLECTIVE )
2208	CALL netcdf_handle_error( 'netcdf_define_header', 448 )
2209	ENDIF
2210	#endif
2211	var_list = TRIM( var_list) // TRIM( do2d(av,i) ) // ';'
2212
2213	ENDIF
2214
2215	i = i + 1
2216
2217	ENDDO
2218
2219	!
2220	!-- No arrays to output. Close the netcdf file and return.
2221	IF ( i == 1 ) RETURN
2222
2223	!
2224	!-- Write the list of variables as global attribute (this is used by
2225	!-- restart runs and by combine_plot_fields)
2226	nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'VAR_LIST', &
2227	var_list )
2228	CALL netcdf_handle_error( 'netcdf_define_header', 121 )
2229
2230	!
2231	!-- Set general no fill, otherwise the performance drops significantly for
2232	!-- parallel output.
2233	nc_stat = NF90_SET_FILL( id_set_xy(av), NF90_NOFILL, oldmode )
2234	CALL netcdf_handle_error( 'netcdf_define_header', 529 )
2235
2236	!
2237	!-- Leave netCDF define mode
2238	nc_stat = NF90_ENDDEF( id_set_xy(av) )
2239	CALL netcdf_handle_error( 'netcdf_define_header', 122 )
2240
2241	!
2242	!-- These data are only written by PE0 for parallel output to increase
2243	!-- the performance.
2244	IF ( myid == 0 .OR. netcdf_data_format < 5 ) THEN
2245
2246	!
2247	!-- Write axis data: z_xy, x, y
2248	ALLOCATE( netcdf_data(1:ns) )
2249
2250	!
2251	!-- Write zu data
2252	DO i = 1, ns
2253	IF( section(i,1) == -1 ) THEN
2254	netcdf_data(i) = -1.0_wp ! section averaged along z
2255	ELSE
2256	netcdf_data(i) = zu( section(i,1) )
2257	ENDIF
2258	ENDDO
2259	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu_xy(av), &
2260	netcdf_data, start = (/ 1 /), &
2261	count = (/ ns /) )
2262	CALL netcdf_handle_error( 'netcdf_define_header', 123 )
2263
2264	!
2265	!-- Write zw data
2266	DO i = 1, ns
2267	IF( section(i,1) == -1 ) THEN
2268	netcdf_data(i) = -1.0_wp ! section averaged along z
2269	ELSE
2270	netcdf_data(i) = zw( section(i,1) )
2271	ENDIF
2272	ENDDO
2273	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zw_xy(av), &
2274	netcdf_data, start = (/ 1 /), &
2275	count = (/ ns /) )
2276	CALL netcdf_handle_error( 'netcdf_define_header', 124 )
2277
2278	!
2279	!-- Write zs data
2280	IF ( land_surface ) THEN
2281	ns_do = 0
2282	DO i = 1, ns
2283	IF( section(i,1) == -1 ) THEN
2284	netcdf_data(i) = 1.0_wp ! section averaged along z
2285	ns_do = ns_do + 1
2286	ELSEIF ( section(i,1) < nzs ) THEN
2287	netcdf_data(i) = - zs( section(i,1) )
2288	ns_do = ns_do + 1
2289	ENDIF
2290	ENDDO
2291
2292	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zs_xy(av), &
2293	netcdf_data(1:ns_do), start = (/ 1 /), &
2294	count = (/ ns_do /) )
2295	CALL netcdf_handle_error( 'netcdf_define_header', 124 )
2296
2297	ENDIF
2298
2299	!
2300	!-- Write gridpoint number data
2301	netcdf_data(1:ns) = section(1:ns,1)
2302	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_ind_z_xy(av), &
2303	netcdf_data, start = (/ 1 /), &
2304	count = (/ ns /) )
2305	CALL netcdf_handle_error( 'netcdf_define_header', 125 )
2306
2307	DEALLOCATE( netcdf_data )
2308
2309	!
2310	!-- Write the cross section height u, t
2311	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu1_xy(av), &
2312	(/ zu(nzb+1) /), start = (/ 1 /), &
2313	count = (/ 1 /) )
2314	CALL netcdf_handle_error( 'netcdf_define_header', 126 )
2315
2316	!
2317	!-- Write data for x (shifted by +dx/2) and xu axis
2318	ALLOCATE( netcdf_data(0:nx) )
2319
2320	DO i = 0, nx
2321	netcdf_data(i) = ( i + 0.5_wp ) * dx
2322	ENDDO
2323
2324	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_x_xy(av), &
2325	netcdf_data, start = (/ 1 /), &
2326	count = (/ nx+1 /) )
2327	CALL netcdf_handle_error( 'netcdf_define_header', 127 )
2328
2329	DO i = 0, nx
2330	netcdf_data(i) = i * dx
2331	ENDDO
2332
2333	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_xu_xy(av), &
2334	netcdf_data, start = (/ 1 /), &
2335	count = (/ nx+1 /) )
2336	CALL netcdf_handle_error( 'netcdf_define_header', 367 )
2337
2338	DEALLOCATE( netcdf_data )
2339
2340	!
2341	!-- Write data for y (shifted by +dy/2) and yv axis
2342	ALLOCATE( netcdf_data(0:ny+1) )
2343
2344	DO i = 0, ny
2345	netcdf_data(i) = ( i + 0.5_wp ) * dy
2346	ENDDO
2347
2348	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_y_xy(av), &
2349	netcdf_data, start = (/ 1 /), &
2350	count = (/ ny+1 /))
2351	CALL netcdf_handle_error( 'netcdf_define_header', 128 )
2352
2353	DO i = 0, ny
2354	netcdf_data(i) = i * dy
2355	ENDDO
2356
2357	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_yv_xy(av), &
2358	netcdf_data, start = (/ 1 /), &
2359	count = (/ ny+1 /))
2360	CALL netcdf_handle_error( 'netcdf_define_header', 368 )
2361
2362	DEALLOCATE( netcdf_data )
2363
2364	ENDIF
2365
2366	!
2367	!-- In case of non-flat topography write height information. Only for
2368	!-- parallel netcdf output.
2369	IF ( TRIM( topography ) /= 'flat' .AND. &
2370	netcdf_data_format > 4 ) THEN
2371
2372	! IF ( nxr == nx .AND. nyn /= ny ) THEN
2373	! nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av), &
2374	! zu_s_inner(nxl:nxr+1,nys:nyn), &
2375	! start = (/ nxl+1, nys+1 /), &
2376	! count = (/ nxr-nxl+2, nyn-nys+1 /) )
2377	! ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN
2378	! nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av), &
2379	! zu_s_inner(nxl:nxr,nys:nyn+1), &
2380	! start = (/ nxl+1, nys+1 /), &
2381	! count = (/ nxr-nxl+1, nyn-nys+2 /) )
2382	! ELSEIF ( nxr == nx .AND. nyn == ny ) THEN
2383	! nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av), &
2384	! zu_s_inner(nxl:nxr+1,nys:nyn+1), &
2385	! start = (/ nxl+1, nys+1 /), &
2386	! count = (/ nxr-nxl+2, nyn-nys+2 /) )
2387	! ELSE
2388	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av), &
2389	zu_s_inner(nxl:nxr,nys:nyn), &
2390	start = (/ nxl+1, nys+1 /), &
2391	count = (/ nxr-nxl+1, nyn-nys+1 /) )
2392	! ENDIF
2393	CALL netcdf_handle_error( 'netcdf_define_header', 427 )
2394
2395	! IF ( nxr == nx .AND. nyn /= ny ) THEN
2396	! nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av), &
2397	! zw_w_inner(nxl:nxr+1,nys:nyn), &
2398	! start = (/ nxl+1, nys+1 /), &
2399	! count = (/ nxr-nxl+2, nyn-nys+1 /) )
2400	! ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN
2401	! nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av), &
2402	! zw_w_inner(nxl:nxr,nys:nyn+1), &
2403	! start = (/ nxl+1, nys+1 /), &
2404	! count = (/ nxr-nxl+1, nyn-nys+2 /) )
2405	! ELSEIF ( nxr == nx .AND. nyn == ny ) THEN
2406	! nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av), &
2407	! zw_w_inner(nxl:nxr+1,nys:nyn+1), &
2408	! start = (/ nxl+1, nys+1 /), &
2409	! count = (/ nxr-nxl+2, nyn-nys+2 /) )
2410	! ELSE
2411	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av), &
2412	zw_w_inner(nxl:nxr,nys:nyn), &
2413	start = (/ nxl+1, nys+1 /), &
2414	count = (/ nxr-nxl+1, nyn-nys+1 /) )
2415	! ENDIF
2416	CALL netcdf_handle_error( 'netcdf_define_header', 428 )
2417
2418	ENDIF
2419
2420	CASE ( 'xy_ext' )
2421
2422	!
2423	!-- Get the list of variables and compare with the actual run.
2424	!-- First var_list_old has to be reset, since GET_ATT does not assign
2425	!-- trailing blanks.
2426	var_list_old = ' '
2427	nc_stat = NF90_GET_ATT( id_set_xy(av), NF90_GLOBAL, 'VAR_LIST', &
2428	var_list_old )
2429	CALL netcdf_handle_error( 'netcdf_define_header', 129 )
2430
2431	var_list = ';'
2432	i = 1
2433	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2434	IF ( INDEX( do2d(av,i), 'xy' ) /= 0 ) THEN
2435	var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
2436	ENDIF
2437	i = i + 1
2438	ENDDO
2439
2440	IF ( av == 0 ) THEN
2441	var = '(xy)'
2442	ELSE
2443	var = '(xy_av)'
2444	ENDIF
2445
2446	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
2447	message_string = 'netCDF file for cross-sections ' // &
2448	TRIM( var ) // ' from previous run found,' // &
2449	'& but this file cannot be extended due to' // &
2450	' variable mismatch.' // &
2451	'&New file is created instead.'
2452	CALL message( 'define_netcdf_header', 'PA0249', 0, 1, 0, 6, 0 )
2453	extend = .FALSE.
2454	RETURN
2455	ENDIF
2456
2457	!
2458	!-- Calculate the number of current sections
2459	ns = 1
2460	DO WHILE ( section(ns,1) /= -9999 .AND. ns <= 100 )
2461	ns = ns + 1
2462	ENDDO
2463	ns = ns - 1
2464
2465	!
2466	!-- Get and compare the number of horizontal cross sections
2467	nc_stat = NF90_INQ_VARID( id_set_xy(av), 'zu_xy', id_var_zu_xy(av) )
2468	CALL netcdf_handle_error( 'netcdf_define_header', 130 )
2469
2470	nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_zu_xy(av), &
2471	dimids = id_dim_zu_xy_old )
2472	CALL netcdf_handle_error( 'netcdf_define_header', 131 )
2473	id_dim_zu_xy(av) = id_dim_zu_xy_old(1)
2474
2475	nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_zu_xy(av), &
2476	len = ns_old )
2477	CALL netcdf_handle_error( 'netcdf_define_header', 132 )
2478
2479	IF ( ns /= ns_old ) THEN
2480	message_string = 'netCDF file for cross-sections ' // &
2481	TRIM( var ) // ' from previous run found,' // &
2482	'&but this file cannot be extended due to' // &
2483	' mismatch in number of' // &
2484	'&cross sections.' // &
2485	'&New file is created instead.'
2486	CALL message( 'define_netcdf_header', 'PA0250', 0, 1, 0, 6, 0 )
2487	extend = .FALSE.
2488	RETURN
2489	ENDIF
2490
2491	!
2492	!-- Get and compare the heights of the cross sections
2493	ALLOCATE( netcdf_data(1:ns_old) )
2494
2495	nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_zu_xy(av), netcdf_data )
2496	CALL netcdf_handle_error( 'netcdf_define_header', 133 )
2497
2498	DO i = 1, ns
2499	IF ( section(i,1) /= -1 ) THEN
2500	IF ( zu(section(i,1)) /= netcdf_data(i) ) THEN
2501	message_string = 'netCDF file for cross-sections ' // &
2502	TRIM( var ) // ' from previous run found,' // &
2503	'&but this file cannot be extended' // &
2504	' due to mismatch in cross' // &
2505	'&section levels.' // &
2506	'&New file is created instead.'
2507	CALL message( 'define_netcdf_header', 'PA0251', &
2508	0, 1, 0, 6, 0 )
2509	extend = .FALSE.
2510	RETURN
2511	ENDIF
2512	ELSE
2513	IF ( -1.0_wp /= netcdf_data(i) ) THEN
2514	message_string = 'netCDF file for cross-sections ' // &
2515	TRIM( var ) // ' from previous run found,' // &
2516	'&but this file cannot be extended' // &
2517	' due to mismatch in cross' // &
2518	'&section levels.' // &
2519	'&New file is created instead.'
2520	CALL message( 'define_netcdf_header', 'PA0251',&
2521	0, 1, 0, 6, 0 )
2522	extend = .FALSE.
2523	RETURN
2524	ENDIF
2525	ENDIF
2526	ENDDO
2527
2528	DEALLOCATE( netcdf_data )
2529
2530	!
2531	!-- Get the id of the time coordinate (unlimited coordinate) and its
2532	!-- last index on the file. The next time level is do2d..count+1.
2533	!-- The current time must be larger than the last output time
2534	!-- on the file.
2535	nc_stat = NF90_INQ_VARID( id_set_xy(av), 'time', id_var_time_xy(av) )
2536	CALL netcdf_handle_error( 'netcdf_define_header', 134 )
2537
2538	nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_time_xy(av), &
2539	dimids = id_dim_time_old )
2540	CALL netcdf_handle_error( 'netcdf_define_header', 135 )
2541	id_dim_time_xy(av) = id_dim_time_old(1)
2542
2543	nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_time_xy(av), &
2544	len = ntime_count )
2545	CALL netcdf_handle_error( 'netcdf_define_header', 136 )
2546
2547	!
2548	!-- For non-parallel output use the last output time level of the netcdf
2549	!-- file because the time dimension is unlimited. In case of parallel
2550	!-- output the variable ntime_count could get the value of 9*10E36 because
2551	!-- the time dimension is limited.
2552	IF ( netcdf_data_format < 5 ) do2d_xy_time_count(av) = ntime_count
2553
2554	nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_time_xy(av), &
2555	last_time_coordinate, &
2556	start = (/ do2d_xy_time_count(av) /), &
2557	count = (/ 1 /) )
2558	CALL netcdf_handle_error( 'netcdf_define_header', 137 )
2559
2560	IF ( last_time_coordinate(1) >= simulated_time ) THEN
2561	message_string = 'netCDF file for cross sections ' // &
2562	TRIM( var ) // ' from previous run found,' // &
2563	'&but this file cannot be extended becaus' // &
2564	'e the current output time' // &
2565	'&is less or equal than the last output t' // &
2566	'ime on this file.' // &
2567	'&New file is created instead.'
2568	CALL message( 'define_netcdf_header', 'PA0252', 0, 1, 0, 6, 0 )
2569	do2d_xy_time_count(av) = 0
2570	extend = .FALSE.
2571	RETURN
2572	ENDIF
2573
2574	IF ( netcdf_data_format > 4 ) THEN
2575	!
2576	!-- Set needed time levels (ntdim) to
2577	!-- saved time levels + to be saved time levels.
2578	IF ( av == 0 ) THEN
2579	ntdim_2d_xy(0) = do2d_xy_time_count(0) + CEILING( &
2580	( end_time - MAX( skip_time_do2d_xy, &
2581	simulated_time_at_begin ) &
2582	) / dt_do2d_xy )
2583	ELSE
2584	ntdim_2d_xy(1) = do2d_xy_time_count(1) + CEILING( &
2585	( end_time - MAX( skip_time_data_output_av, &
2586	simulated_time_at_begin ) &
2587	) / dt_data_output_av )
2588	ENDIF
2589
2590	!
2591	!-- Check if the needed number of output time levels is increased
2592	!-- compared to the number of time levels in the existing file.
2593	IF ( ntdim_2d_xy(av) > ntime_count ) THEN
2594	message_string = 'netCDF file for cross sections ' // &
2595	TRIM( var ) // ' from previous run found,' // &
2596	'&but this file cannot be extended becaus' // &
2597	'e the number of output time levels&has b' // &
2598	'een increased compared to the previous s' // &
2599	'imulation.' // &
2600	'&New file is created instead.'
2601	CALL message( 'define_netcdf_header', 'PA0389', 0, 1, 0, 6, 0 )
2602	do2d_xy_time_count(av) = 0
2603	extend = .FALSE.
2604	!
2605	!-- Recalculate the needed time levels for the new file.
2606	IF ( av == 0 ) THEN
2607	ntdim_2d_xy(0) = CEILING( &
2608	( end_time - MAX( skip_time_do2d_xy, &
2609	simulated_time_at_begin ) &
2610	) / dt_do2d_xy )
2611	ELSE
2612	ntdim_2d_xy(1) = CEILING( &
2613	( end_time - MAX( skip_time_data_output_av, &
2614	simulated_time_at_begin ) &
2615	) / dt_data_output_av )
2616	ENDIF
2617	RETURN
2618	ENDIF
2619	ENDIF
2620
2621	!
2622	!-- Dataset seems to be extendable.
2623	!-- Now get the variable ids.
2624	i = 1
2625	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2626	IF ( INDEX( do2d(av,i), 'xy' ) /= 0 ) THEN
2627	nc_stat = NF90_INQ_VARID( id_set_xy(av), do2d(av,i), &
2628	id_var_do2d(av,i) )
2629	CALL netcdf_handle_error( 'netcdf_define_header', 138 )
2630	#if defined( __netcdf4_parallel )
2631	!
2632	!-- Set collective io operations for parallel io
2633	IF ( netcdf_data_format > 4 ) THEN
2634	nc_stat = NF90_VAR_PAR_ACCESS( id_set_xy(av), &
2635	id_var_do2d(av,i), &
2636	NF90_COLLECTIVE )
2637	CALL netcdf_handle_error( 'netcdf_define_header', 454 )
2638	ENDIF
2639	#endif
2640	ENDIF
2641	i = i + 1
2642	ENDDO
2643
2644	!
2645	!-- Update the title attribute on file
2646	!-- In order to avoid 'data mode' errors if updated attributes are larger
2647	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
2648	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
2649	!-- performance loss due to data copying; an alternative strategy would be
2650	!-- to ensure equal attribute size in a job chain. Maybe revise later.
2651	IF ( av == 0 ) THEN
2652	time_average_text = ' '
2653	ELSE
2654	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
2655	averaging_interval
2656	ENDIF
2657	nc_stat = NF90_REDEF( id_set_xy(av) )
2658	CALL netcdf_handle_error( 'netcdf_define_header', 431 )
2659	nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'title', &
2660	TRIM( run_description_header ) // &
2661	TRIM( time_average_text ) )
2662	CALL netcdf_handle_error( 'netcdf_define_header', 139 )
2663	nc_stat = NF90_ENDDEF( id_set_xy(av) )
2664	CALL netcdf_handle_error( 'netcdf_define_header', 432 )
2665	message_string = 'netCDF file for cross-sections ' // &
2666	TRIM( var ) // ' from previous run found.' // &
2667	'&This file will be extended.'
2668	CALL message( 'define_netcdf_header', 'PA0253', 0, 0, 0, 6, 0 )
2669
2670
2671	CASE ( 'xz_new' )
2672
2673	!
2674	!-- Define some global attributes of the dataset
2675	nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'Conventions', &
2676	'COARDS' )
2677	CALL netcdf_handle_error( 'netcdf_define_header', 140 )
2678
2679	IF ( av == 0 ) THEN
2680	time_average_text = ' '
2681	ELSE
2682	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
2683	averaging_interval
2684	ENDIF
2685	nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'title', &
2686	TRIM( run_description_header ) // &
2687	TRIM( time_average_text ) )
2688	CALL netcdf_handle_error( 'netcdf_define_header', 141 )
2689	IF ( av == 1 ) THEN
2690	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval
2691	nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'time_avg', &
2692	TRIM( time_average_text ) )
2693	CALL netcdf_handle_error( 'netcdf_define_header', 141 )
2694	ENDIF
2695
2696	!
2697	!-- Define time coordinate for xz sections.
2698	!-- For parallel output the time dimensions has to be limited, otherwise
2699	!-- the performance drops significantly.
2700	IF ( netcdf_data_format < 5 ) THEN
2701	CALL netcdf_create_dim( id_set_xz(av), 'time', NF90_UNLIMITED, &
2702	id_dim_time_xz(av), 142 )
2703	ELSE
2704	CALL netcdf_create_dim( id_set_xz(av), 'time', ntdim_2d_xz(av), &
2705	id_dim_time_xz(av), 525 )
2706	ENDIF
2707
2708	CALL netcdf_create_var( id_set_xz(av), (/ id_dim_time_xz(av) /), &
2709	'time', NF90_DOUBLE, id_var_time_xz(av), &
2710	'seconds', '', 143, 144, 000 )
2711	!
2712	!-- Define the spatial dimensions and coordinates for xz-sections.
2713	!-- First, determine the number of vertical sections to be written.
2714	IF ( section(1,2) == -9999 ) THEN
2715	RETURN
2716	ELSE
2717	ns = 1
2718	DO WHILE ( section(ns,2) /= -9999 .AND. ns <= 100 )
2719	ns = ns + 1
2720	ENDDO
2721	ns = ns - 1
2722	ENDIF
2723
2724	!
2725	!-- Define y-axis (for scalar position)
2726	CALL netcdf_create_dim( id_set_xz(av), 'y_xz', ns, id_dim_y_xz(av), &
2727	145 )
2728	CALL netcdf_create_var( id_set_xz(av), (/ id_dim_y_xz(av) /), &
2729	'y_xz', NF90_DOUBLE, id_var_y_xz(av), &
2730	'meters', '', 146, 147, 000 )
2731	!
2732	!-- Define y-axis (for v position)
2733	CALL netcdf_create_dim( id_set_xz(av), 'yv_xz', ns, &
2734	id_dim_yv_xz(av), 369 )
2735	CALL netcdf_create_var( id_set_xz(av), (/ id_dim_yv_xz(av) /), &
2736	'yv_xz', NF90_DOUBLE, id_var_yv_xz(av), &
2737	'meters', '', 370, 371, 000 )
2738	!
2739	!-- Define a variable to store the layer indices of the vertical cross
2740	!-- sections
2741	CALL netcdf_create_var( id_set_xz(av), (/ id_dim_y_xz(av) /), &
2742	'ind_y_xz', NF90_DOUBLE, &
2743	id_var_ind_y_xz(av), 'gridpoints', '', 148, &
2744	149, 000 )
2745	!
2746	!-- Define x-axis (for scalar position)
2747	CALL netcdf_create_dim( id_set_xz(av), 'x', nx+1, id_dim_x_xz(av), &
2748	150 )
2749	CALL netcdf_create_var( id_set_xz(av), (/ id_dim_x_xz(av) /), 'x', &
2750	NF90_DOUBLE, id_var_x_xz(av), 'meters', '', &
2751	151, 152, 000 )
2752	!
2753	!-- Define x-axis (for u position)
2754	CALL netcdf_create_dim( id_set_xz(av), 'xu', nx+1, id_dim_xu_xz(av), &
2755	372 )
2756	CALL netcdf_create_var( id_set_xz(av), (/ id_dim_xu_xz(av) /), 'xu', &
2757	NF90_DOUBLE, id_var_xu_xz(av), 'meters', '', &
2758	373, 374, 000 )
2759	!
2760	!-- Define the three z-axes (zu, zw, and zs)
2761	CALL netcdf_create_dim( id_set_xz(av), 'zu', nz+2, id_dim_zu_xz(av), &
2762	153 )
2763	CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zu_xz(av) /), 'zu', &
2764	NF90_DOUBLE, id_var_zu_xz(av), 'meters', '', &
2765	154, 155, 000 )
2766	CALL netcdf_create_dim( id_set_xz(av), 'zw', nz+2, id_dim_zw_xz(av), &
2767	156 )
2768	CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zw_xz(av) /), 'zw', &
2769	NF90_DOUBLE, id_var_zw_xz(av), 'meters', '', &
2770	157, 158, 000 )
2771
2772	IF ( land_surface ) THEN
2773
2774	CALL netcdf_create_dim( id_set_xz(av), 'zs', nzs, &
2775	id_dim_zs_xz(av), 542 )
2776	CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zs_xz(av) /), &
2777	'zs', NF90_DOUBLE, id_var_zs_xz(av), &
2778	'meters', '', 543, 544, 000 )
2779
2780	ENDIF
2781
2782	!
2783	!-- Define the variables
2784	var_list = ';'
2785	i = 1
2786
2787	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2788
2789	IF ( INDEX( do2d(av,i), 'xz' ) /= 0 ) THEN
2790
2791	!
2792	!-- Check for the grid
2793	found = .FALSE.
2794	SELECT CASE ( do2d(av,i) )
2795	!
2796	!-- Most variables are defined on the zu grid
2797	CASE ( 'e_xz', 'lpt_xz', 'nc_xz', 'nr_xz', 'p_xz', 'pc_xz', &
2798	'pr_xz', 'prr_xz', 'pt_xz', 'q_xz', 'qc_xz', &
2799	'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz', &
2800	'qv_xz', 'rho_ocean_xz', 's_xz', 'sa_xz', 'vpt_xz' )
2801
2802	grid_x = 'x'
2803	grid_y = 'y'
2804	grid_z = 'zu'
2805	!
2806	!-- u grid
2807	CASE ( 'u_xz' )
2808
2809	grid_x = 'xu'
2810	grid_y = 'y'
2811	grid_z = 'zu'
2812	!
2813	!-- v grid
2814	CASE ( 'v_xz' )
2815
2816	grid_x = 'x'
2817	grid_y = 'yv'
2818	grid_z = 'zu'
2819	!
2820	!-- w grid
2821	CASE ( 'w_xz' )
2822
2823	grid_x = 'x'
2824	grid_y = 'y'
2825	grid_z = 'zw'
2826
2827	CASE DEFAULT
2828
2829	!
2830	!-- Check for land surface quantities
2831	IF ( land_surface ) THEN
2832	CALL lsm_define_netcdf_grid( do2d(av,i), found, &
2833	grid_x, grid_y, grid_z )
2834	ENDIF
2835
2836	IF ( .NOT. found ) THEN
2837	CALL tcm_define_netcdf_grid( do2d(av,i), found, &
2838	grid_x, grid_y, grid_z )
2839	ENDIF
2840
2841	!
2842	!-- Check for radiation quantities
2843	IF ( .NOT. found .AND. radiation ) THEN
2844	CALL radiation_define_netcdf_grid( do2d(av,i), found, &
2845	grid_x, grid_y, &
2846	grid_z )
2847	ENDIF
2848
2849	!
2850	!-- Check for chemistry quantities
2851	#if defined( __chem )
2852	IF ( .NOT. found .AND. air_chemistry ) THEN
2853	CALL chem_define_netcdf_grid( do2d(av,i), found, &
2854	grid_x, grid_y, &
2855	grid_z )
2856	ENDIF
2857	#endif
2858
2859	!
2860	!-- Check for user-defined quantities
2861	IF ( .NOT. found ) THEN
2862	CALL user_define_netcdf_grid( do2d(av,i), found, &
2863	grid_x, grid_y, grid_z )
2864	ENDIF
2865
2866	IF ( .NOT. found ) THEN
2867	WRITE ( message_string, * ) 'no grid defined for', &
2868	' variable ', TRIM( do2d(av,i) )
2869	CALL message( 'define_netcdf_header', 'PA0244', &
2870	0, 1, 0, 6, 0 )
2871	ENDIF
2872
2873	END SELECT
2874
2875	!
2876	!-- Select the respective dimension ids
2877	IF ( grid_x == 'x' ) THEN
2878	id_x = id_dim_x_xz(av)
2879	ELSEIF ( grid_x == 'xu' ) THEN
2880	id_x = id_dim_xu_xz(av)
2881	ENDIF
2882
2883	IF ( grid_y == 'y' ) THEN
2884	id_y = id_dim_y_xz(av)
2885	ELSEIF ( grid_y == 'yv' ) THEN
2886	id_y = id_dim_yv_xz(av)
2887	ENDIF
2888
2889	IF ( grid_z == 'zu' ) THEN
2890	id_z = id_dim_zu_xz(av)
2891	ELSEIF ( grid_z == 'zw' ) THEN
2892	id_z = id_dim_zw_xz(av)
2893	ELSEIF ( grid_z == 'zs' ) THEN
2894	id_z = id_dim_zs_xz(av)
2895	ENDIF
2896
2897	!
2898	!-- Define the grid
2899	CALL netcdf_create_var( id_set_xz(av), (/ id_x, id_y, id_z, &
2900	id_dim_time_xz(av) /), do2d(av,i), &
2901	nc_precision(2), id_var_do2d(av,i), &
2902	TRIM( do2d_unit(av,i) ), do2d(av,i), &
2903	159, 160, 355, .TRUE. )
2904
2905	#if defined( __netcdf4_parallel )
2906
2907	IF ( netcdf_data_format > 4 ) THEN
2908	!
2909	!-- Set no fill for every variable to increase performance.
2910	nc_stat = NF90_DEF_VAR_FILL( id_set_xz(av), &
2911	id_var_do2d(av,i), &
2912	1, 0 )
2913	CALL netcdf_handle_error( 'netcdf_define_header', 534 )
2914	!
2915	!-- Set independent io operations for parallel io. Collective io
2916	!-- is only allowed in case of a 1d-decomposition along x,
2917	!-- because otherwise, not all PEs have output data.
2918	IF ( npey == 1 ) THEN
2919	nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av), &
2920	id_var_do2d(av,i), &
2921	NF90_COLLECTIVE )
2922	ELSE
2923	!
2924	!-- Test simulations showed that the output of cross sections
2925	!-- by all PEs in data_output_2d using NF90_COLLECTIVE is
2926	!-- faster than the output by the first row of PEs in
2927	!-- x-direction using NF90_INDEPENDENT.
2928	nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av), &
2929	id_var_do2d(av,i), &
2930	NF90_COLLECTIVE )
2931	! nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av), &
2932	! id_var_do2d(av,i), &
2933	! NF90_INDEPENDENT )
2934	ENDIF
2935	CALL netcdf_handle_error( 'netcdf_define_header', 449 )
2936	ENDIF
2937	#endif
2938	var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
2939
2940	ENDIF
2941
2942	i = i + 1
2943
2944	ENDDO
2945
2946	!
2947	!-- No arrays to output. Close the netcdf file and return.
2948	IF ( i == 1 ) RETURN
2949
2950	!
2951	!-- Write the list of variables as global attribute (this is used by
2952	!-- restart runs and by combine_plot_fields)
2953	nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'VAR_LIST', &
2954	var_list )
2955	CALL netcdf_handle_error( 'netcdf_define_header', 161 )
2956
2957	!
2958	!-- Set general no fill, otherwise the performance drops significantly for
2959	!-- parallel output.
2960	nc_stat = NF90_SET_FILL( id_set_xz(av), NF90_NOFILL, oldmode )
2961	CALL netcdf_handle_error( 'netcdf_define_header', 530 )
2962
2963	!
2964	!-- Leave netCDF define mode
2965	nc_stat = NF90_ENDDEF( id_set_xz(av) )
2966	CALL netcdf_handle_error( 'netcdf_define_header', 162 )
2967
2968	!
2969	!-- These data are only written by PE0 for parallel output to increase
2970	!-- the performance.
2971	IF ( myid == 0 .OR. netcdf_data_format < 5 ) THEN
2972
2973	!
2974	!-- Write axis data: y_xz, x, zu, zw
2975	ALLOCATE( netcdf_data(1:ns) )
2976
2977	!
2978	!-- Write y_xz data (shifted by +dy/2)
2979	DO i = 1, ns
2980	IF( section(i,2) == -1 ) THEN
2981	netcdf_data(i) = -1.0_wp ! section averaged along y
2982	ELSE
2983	netcdf_data(i) = ( section(i,2) + 0.5_wp ) * dy
2984	ENDIF
2985	ENDDO
2986	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_y_xz(av), &
2987	netcdf_data, start = (/ 1 /), &
2988	count = (/ ns /) )
2989	CALL netcdf_handle_error( 'netcdf_define_header', 163 )
2990
2991	!
2992	!-- Write yv_xz data
2993	DO i = 1, ns
2994	IF( section(i,2) == -1 ) THEN
2995	netcdf_data(i) = -1.0_wp ! section averaged along y
2996	ELSE
2997	netcdf_data(i) = section(i,2) * dy
2998	ENDIF
2999	ENDDO
3000	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_yv_xz(av), &
3001	netcdf_data, start = (/ 1 /), &
3002	count = (/ ns /) )
3003	CALL netcdf_handle_error( 'netcdf_define_header', 375 )
3004
3005	!
3006	!-- Write gridpoint number data
3007	netcdf_data(1:ns) = section(1:ns,2)
3008	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_ind_y_xz(av), &
3009	netcdf_data, start = (/ 1 /), &
3010	count = (/ ns /) )
3011	CALL netcdf_handle_error( 'netcdf_define_header', 164 )
3012
3013
3014	DEALLOCATE( netcdf_data )
3015
3016	!
3017	!-- Write data for x (shifted by +dx/2) and xu axis
3018	ALLOCATE( netcdf_data(0:nx) )
3019
3020	DO i = 0, nx
3021	netcdf_data(i) = ( i + 0.5_wp ) * dx
3022	ENDDO
3023
3024	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_x_xz(av), &
3025	netcdf_data, start = (/ 1 /), &
3026	count = (/ nx+1 /) )
3027	CALL netcdf_handle_error( 'netcdf_define_header', 165 )
3028
3029	DO i = 0, nx
3030	netcdf_data(i) = i * dx
3031	ENDDO
3032
3033	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_xu_xz(av), &
3034	netcdf_data, start = (/ 1 /), &
3035	count = (/ nx+1 /) )
3036	CALL netcdf_handle_error( 'netcdf_define_header', 377 )
3037
3038	DEALLOCATE( netcdf_data )
3039
3040	!
3041	!-- Write zu and zw data (vertical axes)
3042	ALLOCATE( netcdf_data(0:nz+1) )
3043
3044	netcdf_data(0:nz+1) = zu(nzb:nzt+1)
3045	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zu_xz(av), &
3046	netcdf_data, start = (/ 1 /), &
3047	count = (/ nz+2 /) )
3048	CALL netcdf_handle_error( 'netcdf_define_header', 166 )
3049
3050	netcdf_data(0:nz+1) = zw(nzb:nzt+1)
3051	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zw_xz(av), &
3052	netcdf_data, start = (/ 1 /), &
3053	count = (/ nz+2 /) )
3054	CALL netcdf_handle_error( 'netcdf_define_header', 167 )
3055
3056	!
3057	!-- Write zs data
3058	IF ( land_surface ) THEN
3059	netcdf_data(0:nzs-1) = - zs(nzb_soil:nzt_soil)
3060	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zs_xz(av), &
3061	netcdf_data(0:nzs), start = (/ 1 /), &
3062	count = (/ nzt_soil-nzb_soil+1 /) )
3063	CALL netcdf_handle_error( 'netcdf_define_header', 548 )
3064	ENDIF
3065
3066
3067	DEALLOCATE( netcdf_data )
3068
3069	ENDIF
3070
3071
3072	CASE ( 'xz_ext' )
3073
3074	!
3075	!-- Get the list of variables and compare with the actual run.
3076	!-- First var_list_old has to be reset, since GET_ATT does not assign
3077	!-- trailing blanks.
3078	var_list_old = ' '
3079	nc_stat = NF90_GET_ATT( id_set_xz(av), NF90_GLOBAL, 'VAR_LIST', &
3080	var_list_old )
3081	CALL netcdf_handle_error( 'netcdf_define_header', 168 )
3082
3083	var_list = ';'
3084	i = 1
3085	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
3086	IF ( INDEX( do2d(av,i), 'xz' ) /= 0 ) THEN
3087	var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
3088	ENDIF
3089	i = i + 1
3090	ENDDO
3091
3092	IF ( av == 0 ) THEN
3093	var = '(xz)'
3094	ELSE
3095	var = '(xz_av)'
3096	ENDIF
3097
3098	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
3099	message_string = 'netCDF file for cross-sections ' // &
3100	TRIM( var ) // ' from previous run found,' // &
3101	'& but this file cannot be extended due to' // &
3102	' variable mismatch.' // &
3103	'&New file is created instead.'
3104	CALL message( 'define_netcdf_header', 'PA0249', 0, 1, 0, 6, 0 )
3105	extend = .FALSE.
3106	RETURN
3107	ENDIF
3108
3109	!
3110	!-- Calculate the number of current sections
3111	ns = 1
3112	DO WHILE ( section(ns,2) /= -9999 .AND. ns <= 100 )
3113	ns = ns + 1
3114	ENDDO
3115	ns = ns - 1
3116
3117	!
3118	!-- Get and compare the number of vertical cross sections
3119	nc_stat = NF90_INQ_VARID( id_set_xz(av), 'y_xz', id_var_y_xz(av) )
3120	CALL netcdf_handle_error( 'netcdf_define_header', 169 )
3121
3122	nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_y_xz(av), &
3123	dimids = id_dim_y_xz_old )
3124	CALL netcdf_handle_error( 'netcdf_define_header', 170 )
3125	id_dim_y_xz(av) = id_dim_y_xz_old(1)
3126
3127	nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_y_xz(av), &
3128	len = ns_old )
3129	CALL netcdf_handle_error( 'netcdf_define_header', 171 )
3130
3131	IF ( ns /= ns_old ) THEN
3132	message_string = 'netCDF file for cross-sections ' // &
3133	TRIM( var ) // ' from previous run found,' // &
3134	'&but this file cannot be extended due to' // &
3135	' mismatch in number of' // &
3136	'&cross sections.' // &
3137	'&New file is created instead.'
3138	CALL message( 'define_netcdf_header', 'PA0250', 0, 1, 0, 6, 0 )
3139	extend = .FALSE.
3140	RETURN
3141	ENDIF
3142
3143	!
3144	!-- Get and compare the heights of the cross sections
3145	ALLOCATE( netcdf_data(1:ns_old) )
3146
3147	nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_y_xz(av), netcdf_data )
3148	CALL netcdf_handle_error( 'netcdf_define_header', 172 )
3149
3150	DO i = 1, ns
3151	IF ( section(i,2) /= -1 ) THEN
3152	IF ( ( ( section(i,2) + 0.5 ) * dy ) /= netcdf_data(i) ) THEN
3153	message_string = 'netCDF file for cross-sections ' // &
3154	TRIM( var ) // ' from previous run found,' // &
3155	'&but this file cannot be extended' // &
3156	' due to mismatch in cross' // &
3157	'&section levels.' // &
3158	'&New file is created instead.'
3159	CALL message( 'define_netcdf_header', 'PA0251', &
3160	0, 1, 0, 6, 0 )
3161	extend = .FALSE.
3162	RETURN
3163	ENDIF
3164	ELSE
3165	IF ( -1.0_wp /= netcdf_data(i) ) THEN
3166	message_string = 'netCDF file for cross-sections ' // &
3167	TRIM( var ) // ' from previous run found,' // &
3168	'&but this file cannot be extended' // &
3169	' due to mismatch in cross' // &
3170	'&section levels.' // &
3171	'&New file is created instead.'
3172	CALL message( 'define_netcdf_header', 'PA0251', &
3173	0, 1, 0, 6, 0 )
3174	extend = .FALSE.
3175	RETURN
3176	ENDIF
3177	ENDIF
3178	ENDDO
3179
3180	DEALLOCATE( netcdf_data )
3181
3182	!
3183	!-- Get the id of the time coordinate (unlimited coordinate) and its
3184	!-- last index on the file. The next time level is do2d..count+1.
3185	!-- The current time must be larger than the last output time
3186	!-- on the file.
3187	nc_stat = NF90_INQ_VARID( id_set_xz(av), 'time', id_var_time_xz(av) )
3188	CALL netcdf_handle_error( 'netcdf_define_header', 173 )
3189
3190	nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_time_xz(av), &
3191	dimids = id_dim_time_old )
3192	CALL netcdf_handle_error( 'netcdf_define_header', 174 )
3193	id_dim_time_xz(av) = id_dim_time_old(1)
3194
3195	nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_time_xz(av), &
3196	len = ntime_count )
3197	CALL netcdf_handle_error( 'netcdf_define_header', 175 )
3198
3199	!
3200	!-- For non-parallel output use the last output time level of the netcdf
3201	!-- file because the time dimension is unlimited. In case of parallel
3202	!-- output the variable ntime_count could get the value of 9*10E36 because
3203	!-- the time dimension is limited.
3204	IF ( netcdf_data_format < 5 ) do2d_xz_time_count(av) = ntime_count
3205
3206	nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_time_xz(av), &
3207	last_time_coordinate, &
3208	start = (/ do2d_xz_time_count(av) /), &
3209	count = (/ 1 /) )
3210	CALL netcdf_handle_error( 'netcdf_define_header', 176 )
3211
3212	IF ( last_time_coordinate(1) >= simulated_time ) THEN
3213	message_string = 'netCDF file for cross sections ' // &
3214	TRIM( var ) // ' from previous run found,' // &
3215	'&but this file cannot be extended becaus' // &
3216	'e the current output time' // &
3217	'&is less or equal than the last output t' // &
3218	'ime on this file.' // &
3219	'&New file is created instead.'
3220	CALL message( 'define_netcdf_header', 'PA0252', 0, 1, 0, 6, 0 )
3221	do2d_xz_time_count(av) = 0
3222	extend = .FALSE.
3223	RETURN
3224	ENDIF
3225
3226	IF ( netcdf_data_format > 4 ) THEN
3227	!
3228	!-- Set needed time levels (ntdim) to
3229	!-- saved time levels + to be saved time levels.
3230	IF ( av == 0 ) THEN
3231	ntdim_2d_xz(0) = do2d_xz_time_count(0) + CEILING( &
3232	( end_time - MAX( skip_time_do2d_xz, &
3233	simulated_time_at_begin ) &
3234	) / dt_do2d_xz )
3235	ELSE
3236	ntdim_2d_xz(1) = do2d_xz_time_count(1) + CEILING( &
3237	( end_time - MAX( skip_time_data_output_av, &
3238	simulated_time_at_begin ) &
3239	) / dt_data_output_av )
3240	ENDIF
3241
3242	!
3243	!-- Check if the needed number of output time levels is increased
3244	!-- compared to the number of time levels in the existing file.
3245	IF ( ntdim_2d_xz(av) > ntime_count ) THEN
3246	message_string = 'netCDF file for cross sections ' // &
3247	TRIM( var ) // ' from previous run found,' // &
3248	'&but this file cannot be extended becaus' // &
3249	'e the number of output time levels&has b' // &
3250	'een increased compared to the previous s' // &
3251	'imulation.' // &
3252	'&New file is created instead.'
3253	CALL message( 'define_netcdf_header', 'PA0390', 0, 1, 0, 6, 0 )
3254	do2d_xz_time_count(av) = 0
3255	extend = .FALSE.
3256	!
3257	!-- Recalculate the needed time levels for the new file.
3258	IF ( av == 0 ) THEN
3259	ntdim_2d_xz(0) = CEILING( &
3260	( end_time - MAX( skip_time_do2d_xz, &
3261	simulated_time_at_begin ) &
3262	) / dt_do2d_xz )
3263	ELSE
3264	ntdim_2d_xz(1) = CEILING( &
3265	( end_time - MAX( skip_time_data_output_av, &
3266	simulated_time_at_begin ) &
3267	) / dt_data_output_av )
3268	ENDIF
3269	RETURN
3270	ENDIF
3271	ENDIF
3272
3273	!
3274	!-- Dataset seems to be extendable.
3275	!-- Now get the variable ids.
3276	i = 1
3277	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
3278	IF ( INDEX( do2d(av,i), 'xz' ) /= 0 ) THEN
3279	nc_stat = NF90_INQ_VARID( id_set_xz(av), do2d(av,i), &
3280	id_var_do2d(av,i) )
3281	CALL netcdf_handle_error( 'netcdf_define_header', 177 )
3282	#if defined( __netcdf4_parallel )
3283	!
3284	!-- Set independent io operations for parallel io. Collective io
3285	!-- is only allowed in case of a 1d-decomposition along x, because
3286	!-- otherwise, not all PEs have output data.
3287	IF ( netcdf_data_format > 4 ) THEN
3288	IF ( npey == 1 ) THEN
3289	nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av), &
3290	id_var_do2d(av,i), &
3291	NF90_COLLECTIVE )
3292	ELSE
3293	!
3294	!-- Test simulations showed that the output of cross sections
3295	!-- by all PEs in data_output_2d using NF90_COLLECTIVE is
3296	!-- faster than the output by the first row of PEs in
3297	!-- x-direction using NF90_INDEPENDENT.
3298	nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av), &
3299	id_var_do2d(av,i), &
3300	NF90_COLLECTIVE )
3301	! nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av), &
3302	! id_var_do2d(av,i), &
3303	! NF90_INDEPENDENT )
3304	ENDIF
3305	CALL netcdf_handle_error( 'netcdf_define_header', 455 )
3306	ENDIF
3307	#endif
3308	ENDIF
3309	i = i + 1
3310	ENDDO
3311
3312	!
3313	!-- Update the title attribute on file
3314	!-- In order to avoid 'data mode' errors if updated attributes are larger
3315	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
3316	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
3317	!-- performance loss due to data copying; an alternative strategy would be
3318	!-- to ensure equal attribute size in a job chain. Maybe revise later.
3319	IF ( av == 0 ) THEN
3320	time_average_text = ' '
3321	ELSE
3322	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
3323	averaging_interval
3324	ENDIF
3325	nc_stat = NF90_REDEF( id_set_xz(av) )
3326	CALL netcdf_handle_error( 'netcdf_define_header', 433 )
3327	nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'title', &
3328	TRIM( run_description_header ) // &
3329	TRIM( time_average_text ) )
3330	CALL netcdf_handle_error( 'netcdf_define_header', 178 )
3331	nc_stat = NF90_ENDDEF( id_set_xz(av) )
3332	CALL netcdf_handle_error( 'netcdf_define_header', 434 )
3333	message_string = 'netCDF file for cross-sections ' // &
3334	TRIM( var ) // ' from previous run found.' // &
3335	'&This file will be extended.'
3336	CALL message( 'define_netcdf_header', 'PA0253', 0, 0, 0, 6, 0 )
3337
3338
3339	CASE ( 'yz_new' )
3340
3341	!
3342	!-- Define some global attributes of the dataset
3343	nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'Conventions', &
3344	'COARDS' )
3345	CALL netcdf_handle_error( 'netcdf_define_header', 179 )
3346
3347	IF ( av == 0 ) THEN
3348	time_average_text = ' '
3349	ELSE
3350	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
3351	averaging_interval
3352	ENDIF
3353	nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'title', &
3354	TRIM( run_description_header ) // &
3355	TRIM( time_average_text ) )
3356	CALL netcdf_handle_error( 'netcdf_define_header', 180 )
3357	IF ( av == 1 ) THEN
3358	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval
3359	nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'time_avg', &
3360	TRIM( time_average_text ) )
3361	CALL netcdf_handle_error( 'netcdf_define_header', 180 )
3362	ENDIF
3363
3364	!
3365	!-- Define time coordinate for yz sections.
3366	!-- For parallel output the time dimensions has to be limited, otherwise
3367	!-- the performance drops significantly.
3368	IF ( netcdf_data_format < 5 ) THEN
3369	CALL netcdf_create_dim( id_set_yz(av), 'time', NF90_UNLIMITED, &
3370	id_dim_time_yz(av), 181 )
3371	ELSE
3372	CALL netcdf_create_dim( id_set_yz(av), 'time', ntdim_2d_yz(av), &
3373	id_dim_time_yz(av), 526 )
3374	ENDIF
3375
3376	CALL netcdf_create_var( id_set_yz(av), (/ id_dim_time_yz(av) /), &
3377	'time', NF90_DOUBLE, id_var_time_yz(av), &
3378	'seconds', '', 182, 183, 000 )
3379	!
3380	!-- Define the spatial dimensions and coordinates for yz-sections.
3381	!-- First, determine the number of vertical sections to be written.
3382	IF ( section(1,3) == -9999 ) THEN
3383	RETURN
3384	ELSE
3385	ns = 1
3386	DO WHILE ( section(ns,3) /= -9999 .AND. ns <= 100 )
3387	ns = ns + 1
3388	ENDDO
3389	ns = ns - 1
3390	ENDIF
3391
3392	!
3393	!-- Define x axis (for scalar position)
3394	CALL netcdf_create_dim( id_set_yz(av), 'x_yz', ns, id_dim_x_yz(av), &
3395	184 )
3396	CALL netcdf_create_var( id_set_yz(av), (/ id_dim_x_yz(av) /), &
3397	'x_yz', NF90_DOUBLE, id_var_x_yz(av), &
3398	'meters', '', 185, 186, 000 )
3399	!
3400	!-- Define x axis (for u position)
3401	CALL netcdf_create_dim( id_set_yz(av), 'xu_yz', ns, &
3402	id_dim_xu_yz(av), 377 )
3403	CALL netcdf_create_var( id_set_yz(av), (/ id_dim_xu_yz(av) /), &
3404	'xu_yz', NF90_DOUBLE, id_var_xu_yz(av), &
3405	'meters', '', 378, 379, 000 )
3406	!
3407	!-- Define a variable to store the layer indices of the vertical cross
3408	!-- sections
3409	CALL netcdf_create_var( id_set_yz(av), (/ id_dim_x_yz(av) /), &
3410	'ind_x_yz', NF90_DOUBLE, &
3411	id_var_ind_x_yz(av), 'gridpoints', '', 187, &
3412	188, 000 )
3413	!
3414	!-- Define y-axis (for scalar position)
3415	CALL netcdf_create_dim( id_set_yz(av), 'y', ny+1, id_dim_y_yz(av), &
3416	189 )
3417	CALL netcdf_create_var( id_set_yz(av), (/ id_dim_y_yz(av) /), 'y', &
3418	NF90_DOUBLE, id_var_y_yz(av), 'meters', '', &
3419	190, 191, 000 )
3420	!
3421	!-- Define y-axis (for v position)
3422	CALL netcdf_create_dim( id_set_yz(av), 'yv', ny+1, id_dim_yv_yz(av), &
3423	380 )
3424	CALL netcdf_create_var( id_set_yz(av), (/ id_dim_yv_yz(av) /), 'yv', &
3425	NF90_DOUBLE, id_var_yv_yz(av), 'meters', '', &
3426	381, 382, 000 )
3427	!
3428	!-- Define the two z-axes (zu and zw)
3429	CALL netcdf_create_dim( id_set_yz(av), 'zu', nz+2, id_dim_zu_yz(av), &
3430	192 )
3431	CALL netcdf_create_var( id_set_yz(av), (/ id_dim_zu_yz(av) /), 'zu', &
3432	NF90_DOUBLE, id_var_zu_yz(av), 'meters', '', &
3433	193, 194, 000 )
3434
3435	CALL netcdf_create_dim( id_set_yz(av), 'zw', nz+2, id_dim_zw_yz(av), &
3436	195 )
3437	CALL netcdf_create_var( id_set_yz(av), (/ id_dim_zw_yz(av) /), 'zw', &
3438	NF90_DOUBLE, id_var_zw_yz(av), 'meters', '', &
3439	196, 197, 000 )
3440
3441	IF ( land_surface ) THEN
3442
3443	CALL netcdf_create_dim( id_set_yz(av), 'zs', nzs, &
3444	id_dim_zs_yz(av), 545 )
3445	CALL netcdf_create_var( id_set_yz(av), (/ id_dim_zs_yz(av) /), &
3446	'zs', NF90_DOUBLE, id_var_zs_yz(av), &
3447	'meters', '', 546, 547, 000 )
3448
3449	ENDIF
3450
3451	!
3452	!-- Define the variables
3453	var_list = ';'
3454	i = 1
3455
3456	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
3457
3458	IF ( INDEX( do2d(av,i), 'yz' ) /= 0 ) THEN
3459
3460	!
3461	!-- Check for the grid
3462	found = .FALSE.
3463	SELECT CASE ( do2d(av,i) )
3464	!
3465	!-- Most variables are defined on the zu grid
3466	CASE ( 'e_yz', 'lpt_yz', 'nc_yz', 'nr_yz', 'p_yz', 'pc_yz', &
3467	'pr_yz','prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz', &
3468	'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz', &
3469	'rho_ocean_yz', 's_yz', 'sa_yz', 'vpt_yz' )
3470
3471	grid_x = 'x'
3472	grid_y = 'y'
3473	grid_z = 'zu'
3474	!
3475	!-- u grid
3476	CASE ( 'u_yz' )
3477
3478	grid_x = 'xu'
3479	grid_y = 'y'
3480	grid_z = 'zu'
3481	!
3482	!-- v grid
3483	CASE ( 'v_yz' )
3484
3485	grid_x = 'x'
3486	grid_y = 'yv'
3487	grid_z = 'zu'
3488	!
3489	!-- w grid
3490	CASE ( 'w_yz' )
3491
3492	grid_x = 'x'
3493	grid_y = 'y'
3494	grid_z = 'zw'
3495
3496
3497	CASE DEFAULT
3498	!
3499	!-- Check for land surface quantities
3500	IF ( land_surface ) THEN
3501	CALL lsm_define_netcdf_grid( do2d(av,i), found, &
3502	grid_x, grid_y, grid_z )
3503	ENDIF
3504
3505	IF ( .NOT. found ) THEN
3506	CALL tcm_define_netcdf_grid( do2d(av,i), found, &
3507	grid_x, grid_y, grid_z )
3508	ENDIF
3509
3510	!
3511	!-- Check for radiation quantities
3512	IF ( .NOT. found .AND. radiation ) THEN
3513	CALL radiation_define_netcdf_grid( do2d(av,i), found, &
3514	grid_x, grid_y, &
3515	grid_z )
3516	ENDIF
3517
3518	!
3519	!-- Check for chemistry quantities
3520	#if defined( __chem )
3521	IF ( .NOT. found .AND. air_chemistry ) THEN
3522	CALL chem_define_netcdf_grid( do2d(av,i), found, &
3523	grid_x, grid_y, &
3524	grid_z )
3525	ENDIF
3526	#endif
3527	!
3528	!-- Check for user-defined quantities
3529	IF ( .NOT. found ) THEN
3530	CALL user_define_netcdf_grid( do2d(av,i), found, &
3531	grid_x, grid_y, grid_z )
3532	ENDIF
3533
3534	IF ( .NOT. found ) THEN
3535	WRITE ( message_string, * ) 'no grid defined for', &
3536	' variable ', TRIM( do2d(av,i) )
3537	CALL message( 'define_netcdf_header', 'PA0244', &
3538	0, 1, 0, 6, 0 )
3539	ENDIF
3540
3541	END SELECT
3542
3543	!
3544	!-- Select the respective dimension ids
3545	IF ( grid_x == 'x' ) THEN
3546	id_x = id_dim_x_yz(av)
3547	ELSEIF ( grid_x == 'xu' ) THEN
3548	id_x = id_dim_xu_yz(av)
3549	ENDIF
3550
3551	IF ( grid_y == 'y' ) THEN
3552	id_y = id_dim_y_yz(av)
3553	ELSEIF ( grid_y == 'yv' ) THEN
3554	id_y = id_dim_yv_yz(av)
3555	ENDIF
3556
3557	IF ( grid_z == 'zu' ) THEN
3558	id_z = id_dim_zu_yz(av)
3559	ELSEIF ( grid_z == 'zw' ) THEN
3560	id_z = id_dim_zw_yz(av)
3561	ELSEIF ( grid_z == 'zs' ) THEN
3562	id_z = id_dim_zs_yz(av)
3563	ENDIF
3564
3565	!
3566	!-- Define the grid
3567	CALL netcdf_create_var( id_set_yz(av), (/ id_x, id_y, id_z, &
3568	id_dim_time_yz(av) /), do2d(av,i), &
3569	nc_precision(3), id_var_do2d(av,i), &
3570	TRIM( do2d_unit(av,i) ), do2d(av,i), &
3571	198, 199, 356, .TRUE. )
3572
3573	#if defined( __netcdf4_parallel )
3574	IF ( netcdf_data_format > 4 ) THEN
3575	!
3576	!-- Set no fill for every variable to increase performance.
3577	nc_stat = NF90_DEF_VAR_FILL( id_set_yz(av), &
3578	id_var_do2d(av,i), &
3579	1, 0 )
3580	CALL netcdf_handle_error( 'netcdf_define_header', 535 )
3581	!
3582	!-- Set independent io operations for parallel io. Collective io
3583	!-- is only allowed in case of a 1d-decomposition along y,
3584	!-- because otherwise, not all PEs have output data.
3585	IF ( npex == 1 ) THEN
3586	nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av), &
3587	id_var_do2d(av,i), &
3588	NF90_COLLECTIVE )
3589	ELSE
3590	!
3591	!-- Test simulations showed that the output of cross sections
3592	!-- by all PEs in data_output_2d using NF90_COLLECTIVE is
3593	!-- faster than the output by the first row of PEs in
3594	!-- y-direction using NF90_INDEPENDENT.
3595	nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av), &
3596	id_var_do2d(av,i), &
3597	NF90_COLLECTIVE )
3598	! nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av), &
3599	! id_var_do2d(av,i), &
3600	! NF90_INDEPENDENT )
3601	ENDIF
3602	CALL netcdf_handle_error( 'netcdf_define_header', 450 )
3603	ENDIF
3604	#endif
3605	var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
3606
3607	ENDIF
3608
3609	i = i + 1
3610
3611	ENDDO
3612
3613	!
3614	!-- No arrays to output. Close the netcdf file and return.
3615	IF ( i == 1 ) RETURN
3616
3617	!
3618	!-- Write the list of variables as global attribute (this is used by
3619	!-- restart runs and by combine_plot_fields)
3620	nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'VAR_LIST', &
3621	var_list )
3622	CALL netcdf_handle_error( 'netcdf_define_header', 200 )
3623
3624	!
3625	!-- Set general no fill, otherwise the performance drops significantly for
3626	!-- parallel output.
3627	nc_stat = NF90_SET_FILL( id_set_yz(av), NF90_NOFILL, oldmode )
3628	CALL netcdf_handle_error( 'netcdf_define_header', 531 )
3629
3630	!
3631	!-- Leave netCDF define mode
3632	nc_stat = NF90_ENDDEF( id_set_yz(av) )
3633	CALL netcdf_handle_error( 'netcdf_define_header', 201 )
3634
3635	!
3636	!-- These data are only written by PE0 for parallel output to increase
3637	!-- the performance.
3638	IF ( myid == 0 .OR. netcdf_data_format < 5 ) THEN
3639
3640	!
3641	!-- Write axis data: x_yz, y, zu, zw
3642	ALLOCATE( netcdf_data(1:ns) )
3643
3644	!
3645	!-- Write x_yz data (shifted by +dx/2)
3646	DO i = 1, ns
3647	IF( section(i,3) == -1 ) THEN
3648	netcdf_data(i) = -1.0_wp ! section averaged along x
3649	ELSE
3650	netcdf_data(i) = ( section(i,3) + 0.5_wp ) * dx
3651	ENDIF
3652	ENDDO
3653	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_x_yz(av), &
3654	netcdf_data, start = (/ 1 /), &
3655	count = (/ ns /) )
3656	CALL netcdf_handle_error( 'netcdf_define_header', 202 )
3657
3658	!
3659	!-- Write x_yz data (xu grid)
3660	DO i = 1, ns
3661	IF( section(i,3) == -1 ) THEN
3662	netcdf_data(i) = -1.0_wp ! section averaged along x
3663	ELSE
3664	netcdf_data(i) = section(i,3) * dx
3665	ENDIF
3666	ENDDO
3667	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_xu_yz(av), &
3668	netcdf_data, start = (/ 1 /), &
3669	count = (/ ns /) )
3670	CALL netcdf_handle_error( 'netcdf_define_header', 383 )
3671
3672	!
3673	!-- Write gridpoint number data
3674	netcdf_data(1:ns) = section(1:ns,3)
3675	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_ind_x_yz(av), &
3676	netcdf_data, start = (/ 1 /), &
3677	count = (/ ns /) )
3678	CALL netcdf_handle_error( 'netcdf_define_header', 203 )
3679
3680	DEALLOCATE( netcdf_data )
3681
3682	!
3683	!-- Write data for y (shifted by +dy/2) and yv axis
3684	ALLOCATE( netcdf_data(0:ny) )
3685
3686	DO j = 0, ny
3687	netcdf_data(j) = ( j + 0.5_wp ) * dy
3688	ENDDO
3689
3690	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_y_yz(av), &
3691	netcdf_data, start = (/ 1 /), &
3692	count = (/ ny+1 /) )
3693	CALL netcdf_handle_error( 'netcdf_define_header', 204 )
3694
3695	DO j = 0, ny
3696	netcdf_data(j) = j * dy
3697	ENDDO
3698
3699	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_yv_yz(av), &
3700	netcdf_data, start = (/ 1 /), &
3701	count = (/ ny+1 /) )
3702	CALL netcdf_handle_error( 'netcdf_define_header', 384 )
3703
3704	DEALLOCATE( netcdf_data )
3705
3706	!
3707	!-- Write zu and zw data (vertical axes)
3708	ALLOCATE( netcdf_data(0:nz+1) )
3709
3710	netcdf_data(0:nz+1) = zu(nzb:nzt+1)
3711	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zu_yz(av), &
3712	netcdf_data, start = (/ 1 /), &
3713	count = (/ nz+2 /) )
3714	CALL netcdf_handle_error( 'netcdf_define_header', 205 )
3715
3716	netcdf_data(0:nz+1) = zw(nzb:nzt+1)
3717	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zw_yz(av), &
3718	netcdf_data, start = (/ 1 /), &
3719	count = (/ nz+2 /) )
3720	CALL netcdf_handle_error( 'netcdf_define_header', 206 )
3721
3722	DEALLOCATE( netcdf_data )
3723
3724	ENDIF
3725
3726
3727	CASE ( 'yz_ext' )
3728
3729	!
3730	!-- Get the list of variables and compare with the actual run.
3731	!-- First var_list_old has to be reset, since GET_ATT does not assign
3732	!-- trailing blanks.
3733	var_list_old = ' '
3734	nc_stat = NF90_GET_ATT( id_set_yz(av), NF90_GLOBAL, 'VAR_LIST', &
3735	var_list_old )
3736	CALL netcdf_handle_error( 'netcdf_define_header', 207 )
3737
3738	var_list = ';'
3739	i = 1
3740	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
3741	IF ( INDEX( do2d(av,i), 'yz' ) /= 0 ) THEN
3742	var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
3743	ENDIF
3744	i = i + 1
3745	ENDDO
3746
3747	IF ( av == 0 ) THEN
3748	var = '(yz)'
3749	ELSE
3750	var = '(yz_av)'
3751	ENDIF
3752
3753	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
3754	message_string = 'netCDF file for cross-sections ' // &
3755	TRIM( var ) // ' from previous run found,' // &
3756	'& but this file cannot be extended due to' // &
3757	' variable mismatch.' // &
3758	'&New file is created instead.'
3759	CALL message( 'define_netcdf_header', 'PA0249', 0, 1, 0, 6, 0 )
3760	extend = .FALSE.
3761	RETURN
3762	ENDIF
3763
3764	!
3765	!-- Calculate the number of current sections
3766	ns = 1
3767	DO WHILE ( section(ns,3) /= -9999 .AND. ns <= 100 )
3768	ns = ns + 1
3769	ENDDO
3770	ns = ns - 1
3771
3772	!
3773	!-- Get and compare the number of vertical cross sections
3774	nc_stat = NF90_INQ_VARID( id_set_yz(av), 'x_yz', id_var_x_yz(av) )
3775	CALL netcdf_handle_error( 'netcdf_define_header', 208 )
3776
3777	nc_stat = NF90_INQUIRE_VARIABLE( id_set_yz(av), id_var_x_yz(av), &
3778	dimids = id_dim_x_yz_old )
3779	CALL netcdf_handle_error( 'netcdf_define_header', 209 )
3780	id_dim_x_yz(av) = id_dim_x_yz_old(1)
3781
3782	nc_stat = NF90_INQUIRE_DIMENSION( id_set_yz(av), id_dim_x_yz(av), &
3783	len = ns_old )
3784	CALL netcdf_handle_error( 'netcdf_define_header', 210 )
3785
3786	IF ( ns /= ns_old ) THEN
3787	message_string = 'netCDF file for cross-sections ' // &
3788	TRIM( var ) // ' from previous run found,' // &
3789	'&but this file cannot be extended due to' // &
3790	' mismatch in number of' // &
3791	'&cross sections.' // &
3792	'&New file is created instead.'
3793	CALL message( 'define_netcdf_header', 'PA0250', 0, 1, 0, 6, 0 )
3794	extend = .FALSE.
3795	RETURN
3796	ENDIF
3797
3798	!
3799	!-- Get and compare the heights of the cross sections
3800	ALLOCATE( netcdf_data(1:ns_old) )
3801
3802	nc_stat = NF90_GET_VAR( id_set_yz(av), id_var_x_yz(av), netcdf_data )
3803	CALL netcdf_handle_error( 'netcdf_define_header', 211 )
3804
3805	DO i = 1, ns
3806	IF ( section(i,3) /= -1 ) THEN
3807	IF ( ( ( section(i,3) + 0.5 ) * dx ) /= netcdf_data(i) ) THEN
3808	message_string = 'netCDF file for cross-sections ' // &
3809	TRIM( var ) // ' from previous run found,' // &
3810	'&but this file cannot be extended' // &
3811	' due to mismatch in cross' // &
3812	'&section levels.' // &
3813	'&New file is created instead.'
3814	CALL message( 'define_netcdf_header', 'PA0251', &
3815	0, 1, 0, 6, 0 )
3816	extend = .FALSE.
3817	RETURN
3818	ENDIF
3819	ELSE
3820	IF ( -1.0_wp /= netcdf_data(i) ) THEN
3821	message_string = 'netCDF file for cross-sections ' // &
3822	TRIM( var ) // ' from previous run found,' // &
3823	'&but this file cannot be extended' // &
3824	' due to mismatch in cross' // &
3825	'&section levels.' // &
3826	'&New file is created instead.'
3827	CALL message( 'define_netcdf_header', 'PA0251', &
3828	0, 1, 0, 6, 0 )
3829	extend = .FALSE.
3830	RETURN
3831	ENDIF
3832	ENDIF
3833	ENDDO
3834
3835	DEALLOCATE( netcdf_data )
3836
3837	!
3838	!-- Get the id of the time coordinate (unlimited coordinate) and its
3839	!-- last index on the file. The next time level is pl2d..count+1.
3840	!-- The current time must be larger than the last output time
3841	!-- on the file.
3842	nc_stat = NF90_INQ_VARID( id_set_yz(av), 'time', id_var_time_yz(av) )
3843	CALL netcdf_handle_error( 'netcdf_define_header', 212 )
3844
3845	nc_stat = NF90_INQUIRE_VARIABLE( id_set_yz(av), id_var_time_yz(av), &
3846	dimids = id_dim_time_old )
3847	CALL netcdf_handle_error( 'netcdf_define_header', 213 )
3848	id_dim_time_yz(av) = id_dim_time_old(1)
3849
3850	nc_stat = NF90_INQUIRE_DIMENSION( id_set_yz(av), id_dim_time_yz(av), &
3851	len = ntime_count )
3852	CALL netcdf_handle_error( 'netcdf_define_header', 214 )
3853
3854	!
3855	!-- For non-parallel output use the last output time level of the netcdf
3856	!-- file because the time dimension is unlimited. In case of parallel
3857	!-- output the variable ntime_count could get the value of 9*10E36 because
3858	!-- the time dimension is limited.
3859	IF ( netcdf_data_format < 5 ) do2d_yz_time_count(av) = ntime_count
3860
3861	nc_stat = NF90_GET_VAR( id_set_yz(av), id_var_time_yz(av), &
3862	last_time_coordinate, &
3863	start = (/ do2d_yz_time_count(av) /), &
3864	count = (/ 1 /) )
3865	CALL netcdf_handle_error( 'netcdf_define_header', 215 )
3866
3867	IF ( last_time_coordinate(1) >= simulated_time ) THEN
3868	message_string = 'netCDF file for cross sections ' // &
3869	TRIM( var ) // ' from previous run found,' // &
3870	'&but this file cannot be extended becaus' // &
3871	'e the current output time' // &
3872	'&is less or equal than the last output t' // &
3873	'ime on this file.' // &
3874	'&New file is created instead.'
3875	CALL message( 'define_netcdf_header', 'PA0252', 0, 1, 0, 6, 0 )
3876	do2d_yz_time_count(av) = 0
3877	extend = .FALSE.
3878	RETURN
3879	ENDIF
3880
3881	IF ( netcdf_data_format > 4 ) THEN
3882	!
3883	!-- Set needed time levels (ntdim) to
3884	!-- saved time levels + to be saved time levels.
3885	IF ( av == 0 ) THEN
3886	ntdim_2d_yz(0) = do2d_yz_time_count(0) + CEILING( &
3887	( end_time - MAX( skip_time_do2d_yz, &
3888	simulated_time_at_begin ) &
3889	) / dt_do2d_yz )
3890	ELSE
3891	ntdim_2d_yz(1) = do2d_yz_time_count(1) + CEILING( &
3892	( end_time - MAX( skip_time_data_output_av, &
3893	simulated_time_at_begin ) &
3894	) / dt_data_output_av )
3895	ENDIF
3896
3897	!
3898	!-- Check if the needed number of output time levels is increased
3899	!-- compared to the number of time levels in the existing file.
3900	IF ( ntdim_2d_yz(av) > ntime_count ) THEN
3901	message_string = 'netCDF file for cross sections ' // &
3902	TRIM( var ) // ' from previous run found,' // &
3903	'&but this file cannot be extended becaus' // &
3904	'e the number of output time levels&has b' // &
3905	'een increased compared to the previous s' // &
3906	'imulation.' // &
3907	'&New file is created instead.'
3908	CALL message( 'define_netcdf_header', 'PA0391', 0, 1, 0, 6, 0 )
3909	do2d_yz_time_count(av) = 0
3910	extend = .FALSE.
3911	!
3912	!-- Recalculate the needed time levels for the new file.
3913	IF ( av == 0 ) THEN
3914	ntdim_2d_yz(0) = CEILING( &
3915	( end_time - MAX( skip_time_do2d_yz, &
3916	simulated_time_at_begin ) &
3917	) / dt_do2d_yz )
3918	ELSE
3919	ntdim_2d_yz(1) = CEILING( &
3920	( end_time - MAX( skip_time_data_output_av, &
3921	simulated_time_at_begin ) &
3922	) / dt_data_output_av )
3923	ENDIF
3924	RETURN
3925	ENDIF
3926	ENDIF
3927
3928	!
3929	!-- Dataset seems to be extendable.
3930	!-- Now get the variable ids.
3931	i = 1
3932	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
3933	IF ( INDEX( do2d(av,i), 'yz' ) /= 0 ) THEN
3934	nc_stat = NF90_INQ_VARID( id_set_yz(av), do2d(av,i), &
3935	id_var_do2d(av,i) )
3936	CALL netcdf_handle_error( 'netcdf_define_header', 216 )
3937	#if defined( __netcdf4_parallel )
3938	!
3939	!-- Set independent io operations for parallel io. Collective io
3940	!-- is only allowed in case of a 1d-decomposition along y, because
3941	!-- otherwise, not all PEs have output data.
3942	IF ( netcdf_data_format > 4 ) THEN
3943	IF ( npex == 1 ) THEN
3944	nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av), &
3945	id_var_do2d(av,i), &
3946	NF90_COLLECTIVE )
3947	ELSE
3948	!
3949	!-- Test simulations showed that the output of cross sections
3950	!-- by all PEs in data_output_2d using NF90_COLLECTIVE is
3951	!-- faster than the output by the first row of PEs in
3952	!-- y-direction using NF90_INDEPENDENT.
3953	nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av), &
3954	id_var_do2d(av,i), &
3955	NF90_COLLECTIVE )
3956	! nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av), &
3957	! id_var_do2d(av,i), &
3958	! NF90_INDEPENDENT )
3959	ENDIF
3960	CALL netcdf_handle_error( 'netcdf_define_header', 450 )
3961	ENDIF
3962	#endif
3963	ENDIF
3964	i = i + 1
3965	ENDDO
3966
3967	!
3968	!-- Update the title attribute on file
3969	!-- In order to avoid 'data mode' errors if updated attributes are larger
3970	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
3971	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
3972	!-- performance loss due to data copying; an alternative strategy would be
3973	!-- to ensure equal attribute size in a job chain. Maybe revise later.
3974	IF ( av == 0 ) THEN
3975	time_average_text = ' '
3976	ELSE
3977	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
3978	averaging_interval
3979	ENDIF
3980	nc_stat = NF90_REDEF( id_set_yz(av) )
3981	CALL netcdf_handle_error( 'netcdf_define_header', 435 )
3982	nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'title', &
3983	TRIM( run_description_header ) // &
3984	TRIM( time_average_text ) )
3985	CALL netcdf_handle_error( 'netcdf_define_header', 217 )
3986	nc_stat = NF90_ENDDEF( id_set_yz(av) )
3987	CALL netcdf_handle_error( 'netcdf_define_header', 436 )
3988	message_string = 'netCDF file for cross-sections ' // &
3989	TRIM( var ) // ' from previous run found.' // &
3990	'&This file will be extended.'
3991	CALL message( 'define_netcdf_header', 'PA0253', 0, 0, 0, 6, 0 )
3992
3993
3994	CASE ( 'pr_new' )
3995
3996	!
3997	!-- Define some global attributes of the dataset
3998	IF ( averaging_interval_pr /= 0.0_wp ) THEN
3999	WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
4000	averaging_interval_pr
4001	nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'title', &
4002	TRIM( run_description_header ) // &
4003	TRIM( time_average_text ) )
4004	CALL netcdf_handle_error( 'netcdf_define_header', 218 )
4005
4006	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval_pr
4007	nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'time_avg', &
4008	TRIM( time_average_text ) )
4009	ELSE
4010	nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'title', &
4011	TRIM( run_description_header ) )
4012	ENDIF
4013	CALL netcdf_handle_error( 'netcdf_define_header', 219 )
4014
4015	!
4016	!-- Write number of columns and rows of coordinate systems to be plotted
4017	!-- on one page to the netcdf header.
4018	!-- This information can be used by palmplot.
4019	nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, &
4020	'no_rows', &
4021	profile_rows )
4022	CALL netcdf_handle_error( 'netcdf_define_header', 519 )
4023
4024	nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, &
4025	'no_columns', &
4026	profile_columns )
4027	CALL netcdf_handle_error( 'netcdf_define_header', 520 )
4028
4029
4030	cross_profiles_adj = ADJUSTL( cross_profiles )
4031	cross_profiles_numb = 999999
4032	cross_profiles_char = ''
4033
4034	!
4035	!-- Each profile defined in cross_profiles is written to an array
4036	!-- (cross_profiles_char). The number of the respective coordinate
4037	!-- system is assigned in a second array (cross_profiles_numb).
4038	k = 1
4039
4040	DO i = 1, crmax
4041
4042	IF ( TRIM( cross_profiles_adj(i) ) == ' ' ) EXIT
4043	delim_old = 0
4044
4045	DO j = 1, crmax
4046	delim = INDEX( cross_profiles_adj(i)(delim_old+1:), ' ' )
4047	IF ( delim == 1 ) EXIT
4048	kk = MIN( crmax, k )
4049	cross_profiles_char(kk) = cross_profiles_adj(i)(delim_old+1: &
4050	delim_old+delim-1)
4051	cross_profiles_numb(kk) = i
4052	k = k + 1
4053	cross_profiles_maxi = i
4054	delim_old = delim_old + delim
4055	ENDDO
4056
4057	ENDDO
4058
4059	cross_profiles_count = MIN( crmax, k-1 )
4060	!
4061	!-- Check if all profiles defined in cross_profiles are defined in
4062	!-- data_output_pr. If not, they will be skipped.
4063	DO i = 1, cross_profiles_count
4064	DO j = 1, dopr_n
4065
4066	IF ( TRIM(cross_profiles_char(i)) == TRIM(data_output_pr(j)) ) &
4067	THEN
4068	EXIT
4069	ENDIF
4070
4071	IF ( j == dopr_n ) THEN
4072	cross_profiles_numb(i) = 999999
4073	ENDIF
4074
4075	ENDDO
4076	ENDDO
4077
4078	DO i = 1, crmax
4079	IF ( cross_profiles_numb(i) == 999999 ) THEN
4080	DO j = i + 1, crmax
4081	IF ( cross_profiles_numb(j) /= 999999 ) THEN
4082	cross_profiles_char(i) = cross_profiles_char(j)
4083	cross_profiles_numb(i) = cross_profiles_numb(j)
4084	cross_profiles_numb(j) = 999999
4085	EXIT
4086	ENDIF
4087	ENDDO
4088	ENDIF
4089	ENDDO
4090
4091	DO i = 1, crmax-1
4092	IF ( cross_profiles_numb(i + 1) == 999999 ) THEN
4093	cross_profiles_count = i
4094	EXIT
4095	ENDIF
4096	ENDDO
4097	!
4098	!-- Check if all profiles defined in data_output_pr are defined in
4099	!-- cross_profiles. If not, they will be added to cross_profiles.
4100	DO i = 1, dopr_n
4101	DO j = 1, cross_profiles_count
4102
4103	IF ( TRIM(cross_profiles_char(j)) == TRIM(data_output_pr(i))) &
4104	THEN
4105	EXIT
4106	ENDIF
4107
4108	IF (( j == cross_profiles_count ) .AND. &
4109	( cross_profiles_count <= crmax - 1)) THEN
4110	cross_profiles_count = cross_profiles_count + 1
4111	cross_profiles_maxi = cross_profiles_maxi + 1
4112	cross_profiles_char(MIN( crmax, cross_profiles_count )) = &
4113	TRIM( data_output_pr(i) )
4114	cross_profiles_numb(MIN( crmax, cross_profiles_count )) = &
4115	cross_profiles_maxi
4116	ENDIF
4117
4118	ENDDO
4119	ENDDO
4120
4121	IF ( cross_profiles_count >= crmax ) THEN
4122	message_string = 'It is not allowed to arrange more than ' &
4123	// '100 profiles with & cross_profiles. Apart'&
4124	// ' from that, all profiles are saved & to ' &
4125	// 'the netCDF file.'
4126	CALL message( 'define_netcdf_header', 'PA0354', 0, 0, 0, 6, 0 )
4127	ENDIF
4128
4129	!
4130	!-- Writing cross_profiles to netcdf header. This information can be
4131	!-- used by palmplot. Each profile is separated by ",", each cross is
4132	!-- separated by ";".
4133	char_cross_profiles = ';'
4134	id_last = 1
4135	cross_profiles_count = MIN( cross_profiles_count, crmax )
4136
4137	DO i = 1, cross_profiles_count
4138
4139	IF ( cross_profiles_numb(i) /= 999999 ) THEN
4140	IF ( TRIM( char_cross_profiles ) == ';' ) THEN
4141	char_cross_profiles = TRIM( char_cross_profiles ) // &
4142	TRIM( cross_profiles_char(i) )
4143	ELSEIF ( id_last == cross_profiles_numb(i) ) THEN
4144	char_cross_profiles = TRIM( char_cross_profiles ) // &
4145	',' // TRIM( cross_profiles_char(i) )
4146	ELSE
4147	char_cross_profiles = TRIM( char_cross_profiles ) // &
4148	';' // TRIM( cross_profiles_char(i) )
4149	ENDIF
4150	id_last = cross_profiles_numb(i)
4151	ENDIF
4152
4153	ENDDO
4154
4155	char_cross_profiles = TRIM( char_cross_profiles ) // ';'
4156
4157	nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'cross_profiles', &
4158	TRIM( char_cross_profiles ) )
4159	CALL netcdf_handle_error( 'netcdf_define_header', 521 )
4160
4161	!
4162	!-- Define time coordinate for profiles (unlimited dimension)
4163	CALL netcdf_create_dim( id_set_pr, 'time', NF90_UNLIMITED, &
4164	id_dim_time_pr, 220 )
4165	CALL netcdf_create_var( id_set_pr, (/ id_dim_time_pr /), 'time', &
4166	NF90_DOUBLE, id_var_time_pr, 'seconds', '', &
4167	221, 222, 000 )
4168	!
4169	!-- Define the variables
4170	var_list = ';'
4171	DO i = 1, dopr_n
4172
4173	IF ( statistic_regions == 0 ) THEN
4174
4175	!
4176	!-- Define the z-axes (each variable gets its own z-axis)
4177	CALL netcdf_create_dim( id_set_pr, &
4178	'z' // TRIM( data_output_pr(i) ), &
4179	nzt+2-nzb, id_dim_z_pr(i,0), 223 )
4180	CALL netcdf_create_var( id_set_pr, (/ id_dim_z_pr(i,0) /), &
4181	'z' // TRIM( data_output_pr(i) ), &
4182	NF90_DOUBLE, id_var_z_pr(i,0), &
4183	'meters', '', 224, 225, 000 )
4184	!
4185	!-- Define the variable
4186	CALL netcdf_create_var( id_set_pr, (/ id_dim_z_pr(i,0), &
4187	id_dim_time_pr /), data_output_pr(i), &
4188	nc_precision(5), id_var_dopr(i,0), &
4189	TRIM( dopr_unit(i) ), &
4190	TRIM( data_output_pr(i) ), 226, 227, &
4191	228 )
4192
4193	var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) // ';'
4194
4195	ELSE
4196	!
4197	!-- If statistic regions are defined, add suffix _SR+#SR to the
4198	!-- names
4199	DO j = 0, statistic_regions
4200	WRITE ( suffix, '(''_'',I2.2)' ) j
4201
4202	!
4203	!-- Define the z-axes (each variable gets it own z-axis)
4204	CALL netcdf_create_dim( id_set_pr, 'z' // &
4205	TRIM(data_output_pr(i)) // suffix, &
4206	nzt+2-nzb, id_dim_z_pr(i,j), 229 )
4207	CALL netcdf_create_var( id_set_pr, (/ id_dim_z_pr(i,j) /), &
4208	'z' // TRIM(data_output_pr(i)) // &
4209	suffix, NF90_DOUBLE, &
4210	id_var_z_pr(i,j), 'meters', '', &
4211	230, 231, 000 )
4212	!
4213	!-- Define the variable
4214	CALL netcdf_create_var( id_set_pr, (/ id_dim_z_pr(i,j), &
4215	id_dim_time_pr /), &
4216	TRIM(data_output_pr(i)) // suffix, &
4217	nc_precision(5), id_var_dopr(i,j), &
4218	TRIM( dopr_unit(i) ), &
4219	TRIM( data_output_pr(i) ) // &
4220	' SR ', 232, 233, 234 )
4221
4222	var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) // &
4223	suffix // ';'
4224
4225	ENDDO
4226
4227	ENDIF
4228
4229	ENDDO
4230
4231	!
4232	!-- Write the list of variables as global attribute (this is used by
4233	!-- restart runs)
4234	nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'VAR_LIST', var_list )
4235	CALL netcdf_handle_error( 'netcdf_define_header', 235 )
4236
4237	!
4238	!-- Define normalization variables (as time series)
4239	DO i = 1, dopr_norm_num
4240
4241	CALL netcdf_create_var( id_set_pr, (/ id_dim_time_pr /), &
4242	'NORM_' // TRIM( dopr_norm_names(i) ), &
4243	nc_precision(5), id_var_norm_dopr(i), &
4244	'', TRIM( dopr_norm_longnames(i) ), 236, &
4245	000, 237 )
4246
4247	ENDDO
4248
4249	!
4250	!-- Leave netCDF define mode
4251	nc_stat = NF90_ENDDEF( id_set_pr )
4252	CALL netcdf_handle_error( 'netcdf_define_header', 238 )
4253
4254	!
4255	!-- Write z-axes data
4256	DO i = 1, dopr_n
4257	DO j = 0, statistic_regions
4258
4259	nc_stat = NF90_PUT_VAR( id_set_pr, id_var_z_pr(i,j), &
4260	hom(nzb:nzt+1,2,dopr_index(i),0), &
4261	start = (/ 1 /), &
4262	count = (/ nzt-nzb+2 /) )
4263	CALL netcdf_handle_error( 'netcdf_define_header', 239 )
4264
4265	ENDDO
4266	ENDDO
4267
4268
4269	CASE ( 'pr_ext' )
4270
4271	!
4272	!-- Get the list of variables and compare with the actual run.
4273	!-- First var_list_old has to be reset, since GET_ATT does not assign
4274	!-- trailing blanks.
4275	var_list_old = ' '
4276	nc_stat = NF90_GET_ATT( id_set_pr, NF90_GLOBAL, 'VAR_LIST', &
4277	var_list_old )
4278	CALL netcdf_handle_error( 'netcdf_define_header', 240 )
4279
4280	var_list = ';'
4281	DO i = 1, dopr_n
4282
4283	IF ( statistic_regions == 0 ) THEN
4284	var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) // ';'
4285	ELSE
4286	DO j = 0, statistic_regions
4287	WRITE ( suffix, '(''_'',I2.2)' ) j
4288	var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) // &
4289	suffix // ';'
4290	ENDDO
4291	ENDIF
4292
4293	ENDDO
4294
4295	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
4296	message_string = 'netCDF file for vertical profiles ' // &
4297	'from previous run found,' // &
4298	'& but this file cannot be extended due to' // &
4299	' variable mismatch.' // &
4300	'&New file is created instead.'
4301	CALL message( 'define_netcdf_header', 'PA0254', 0, 1, 0, 6, 0 )
4302	extend = .FALSE.
4303	RETURN
4304	ENDIF
4305
4306	!
4307	!-- Get the id of the time coordinate (unlimited coordinate) and its
4308	!-- last index on the file. The next time level is dopr..count+1.
4309	!-- The current time must be larger than the last output time
4310	!-- on the file.
4311	nc_stat = NF90_INQ_VARID( id_set_pr, 'time', id_var_time_pr )
4312	CALL netcdf_handle_error( 'netcdf_define_header', 241 )
4313
4314	nc_stat = NF90_INQUIRE_VARIABLE( id_set_pr, id_var_time_pr, &
4315	dimids = id_dim_time_old )
4316	CALL netcdf_handle_error( 'netcdf_define_header', 242 )
4317	id_dim_time_pr = id_dim_time_old(1)
4318
4319	nc_stat = NF90_INQUIRE_DIMENSION( id_set_pr, id_dim_time_pr, &
4320	len = dopr_time_count )
4321	CALL netcdf_handle_error( 'netcdf_define_header', 243 )
4322
4323	nc_stat = NF90_GET_VAR( id_set_pr, id_var_time_pr, &
4324	last_time_coordinate, &
4325	start = (/ dopr_time_count /), &
4326	count = (/ 1 /) )
4327	CALL netcdf_handle_error( 'netcdf_define_header', 244 )
4328
4329	IF ( last_time_coordinate(1) >= simulated_time ) THEN
4330	message_string = 'netCDF file for vertical profiles ' // &
4331	'from previous run found,' // &
4332	'&but this file cannot be extended becaus' // &
4333	'e the current output time' // &
4334	'&is less or equal than the last output t' // &
4335	'ime on this file.' // &
4336	'&New file is created instead.'
4337	CALL message( 'define_netcdf_header', 'PA0255', 0, 1, 0, 6, 0 )
4338	dopr_time_count = 0
4339	extend = .FALSE.
4340	RETURN
4341	ENDIF
4342
4343	!
4344	!-- Dataset seems to be extendable.
4345	!-- Now get the variable ids.
4346	i = 1
4347	DO i = 1, dopr_n
4348
4349	IF ( statistic_regions == 0 ) THEN
4350	nc_stat = NF90_INQ_VARID( id_set_pr, data_output_pr(i), &
4351	id_var_dopr(i,0) )
4352	CALL netcdf_handle_error( 'netcdf_define_header', 245 )
4353	ELSE
4354	DO j = 0, statistic_regions
4355	WRITE ( suffix, '(''_'',I2.2)' ) j
4356	netcdf_var_name = TRIM( data_output_pr(i) ) // suffix
4357	nc_stat = NF90_INQ_VARID( id_set_pr, netcdf_var_name, &
4358	id_var_dopr(i,j) )
4359	CALL netcdf_handle_error( 'netcdf_define_header', 246 )
4360	ENDDO
4361	ENDIF
4362
4363	ENDDO
4364
4365	!
4366	!-- Get ids of the normalization variables
4367	DO i = 1, dopr_norm_num
4368	nc_stat = NF90_INQ_VARID( id_set_pr, &
4369	'NORM_' // TRIM( dopr_norm_names(i) ), &
4370	id_var_norm_dopr(i) )
4371	CALL netcdf_handle_error( 'netcdf_define_header', 247 )
4372	ENDDO
4373
4374	!
4375	!-- Update the title attribute on file
4376	!-- In order to avoid 'data mode' errors if updated attributes are larger
4377	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
4378	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
4379	!-- performance loss due to data copying; an alternative strategy would be
4380	!-- to ensure equal attribute size in a job chain. Maybe revise later.
4381	IF ( averaging_interval_pr == 0.0_wp ) THEN
4382	time_average_text = ' '
4383	ELSE
4384	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
4385	averaging_interval_pr
4386	ENDIF
4387	nc_stat = NF90_REDEF( id_set_pr )
4388	CALL netcdf_handle_error( 'netcdf_define_header', 437 )
4389	nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'title', &
4390	TRIM( run_description_header ) // &
4391	TRIM( time_average_text ) )
4392	CALL netcdf_handle_error( 'netcdf_define_header', 248 )
4393
4394	nc_stat = NF90_ENDDEF( id_set_pr )
4395	CALL netcdf_handle_error( 'netcdf_define_header', 438 )
4396	message_string = 'netCDF file for vertical profiles ' // &
4397	'from previous run found.' // &
4398	'&This file will be extended.'
4399	CALL message( 'define_netcdf_header', 'PA0256', 0, 0, 0, 6, 0 )
4400
4401
4402	CASE ( 'ts_new' )
4403
4404	!
4405	!-- Define some global attributes of the dataset
4406	nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'title', &
4407	TRIM( run_description_header ) )
4408	CALL netcdf_handle_error( 'netcdf_define_header', 249 )
4409
4410	!
4411	!-- Define time coordinate for time series (unlimited dimension)
4412	CALL netcdf_create_dim( id_set_ts, 'time', NF90_UNLIMITED, &
4413	id_dim_time_ts, 250 )
4414	CALL netcdf_create_var( id_set_ts, (/ id_dim_time_ts /), 'time', &
4415	NF90_DOUBLE, id_var_time_ts, 'seconds', '', &
4416	251, 252, 000 )
4417	!
4418	!-- Define the variables
4419	var_list = ';'
4420	DO i = 1, dots_num
4421
4422	IF ( statistic_regions == 0 ) THEN
4423
4424	CALL netcdf_create_var( id_set_ts, (/ id_dim_time_ts /), &
4425	dots_label(i), nc_precision(6), &
4426	id_var_dots(i,0), &
4427	TRIM( dots_unit(i) ), &
4428	TRIM( dots_label(i) ), 253, 254, 255 )
4429
4430	var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // ';'
4431
4432	ELSE
4433	!
4434	!-- If statistic regions are defined, add suffix _SR+#SR to the
4435	!-- names
4436	DO j = 0, statistic_regions
4437	WRITE ( suffix, '(''_'',I2.2)' ) j
4438
4439	CALL netcdf_create_var( id_set_ts, (/ id_dim_time_ts /), &
4440	TRIM( dots_label(i) ) // suffix, &
4441	nc_precision(6), id_var_dots(i,j), &
4442	TRIM( dots_unit(i) ), &
4443	TRIM( dots_label(i) ) // ' SR ' // &
4444	suffix(2:2), 256, 257, 347)
4445
4446	var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // &
4447	suffix // ';'
4448
4449	ENDDO
4450
4451	ENDIF
4452
4453	ENDDO
4454
4455	!
4456	!-- Write the list of variables as global attribute (this is used by
4457	!-- restart runs)
4458	nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'VAR_LIST', var_list )
4459	CALL netcdf_handle_error( 'netcdf_define_header', 258 )
4460
4461	!
4462	!-- Leave netCDF define mode
4463	nc_stat = NF90_ENDDEF( id_set_ts )
4464	CALL netcdf_handle_error( 'netcdf_define_header', 259 )
4465
4466
4467	CASE ( 'ts_ext' )
4468
4469	!
4470	!-- Get the list of variables and compare with the actual run.
4471	!-- First var_list_old has to be reset, since GET_ATT does not assign
4472	!-- trailing blanks.
4473	var_list_old = ' '
4474	nc_stat = NF90_GET_ATT( id_set_ts, NF90_GLOBAL, 'VAR_LIST', &
4475	var_list_old )
4476	CALL netcdf_handle_error( 'netcdf_define_header', 260 )
4477
4478	var_list = ';'
4479	i = 1
4480	DO i = 1, dots_num
4481
4482	IF ( statistic_regions == 0 ) THEN
4483	var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // ';'
4484	ELSE
4485	DO j = 0, statistic_regions
4486	WRITE ( suffix, '(''_'',I2.2)' ) j
4487	var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // &
4488	suffix // ';'
4489	ENDDO
4490	ENDIF
4491
4492	ENDDO
4493
4494	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
4495	message_string = 'netCDF file for time series ' // &
4496	'from previous run found,' // &
4497	'& but this file cannot be extended due to' // &
4498	' variable mismatch.' // &
4499	'&New file is created instead.'
4500	CALL message( 'define_netcdf_header', 'PA0257', 0, 1, 0, 6, 0 )
4501	extend = .FALSE.
4502	RETURN
4503	ENDIF
4504
4505	!
4506	!-- Get the id of the time coordinate (unlimited coordinate) and its
4507	!-- last index on the file. The next time level is dots..count+1.
4508	!-- The current time must be larger than the last output time
4509	!-- on the file.
4510	nc_stat = NF90_INQ_VARID( id_set_ts, 'time', id_var_time_ts )
4511	CALL netcdf_handle_error( 'netcdf_define_header', 261 )
4512
4513	nc_stat = NF90_INQUIRE_VARIABLE( id_set_ts, id_var_time_ts, &
4514	dimids = id_dim_time_old )
4515	CALL netcdf_handle_error( 'netcdf_define_header', 262 )
4516	id_dim_time_ts = id_dim_time_old(1)
4517
4518	nc_stat = NF90_INQUIRE_DIMENSION( id_set_ts, id_dim_time_ts, &
4519	len = dots_time_count )
4520	CALL netcdf_handle_error( 'netcdf_define_header', 263 )
4521
4522	nc_stat = NF90_GET_VAR( id_set_ts, id_var_time_ts, &
4523	last_time_coordinate, &
4524	start = (/ dots_time_count /), &
4525	count = (/ 1 /) )
4526	CALL netcdf_handle_error( 'netcdf_define_header', 264 )
4527
4528	IF ( last_time_coordinate(1) >= simulated_time ) THEN
4529	message_string = 'netCDF file for time series ' // &
4530	'from previous run found,' // &
4531	'&but this file cannot be extended becaus' // &
4532	'e the current output time' // &
4533	'&is less or equal than the last output t' // &
4534	'ime on this file.' // &
4535	'&New file is created instead.'
4536	CALL message( 'define_netcdf_header', 'PA0258', 0, 1, 0, 6, 0 )
4537	dots_time_count = 0
4538	extend = .FALSE.
4539	RETURN
4540	ENDIF
4541
4542	!
4543	!-- Dataset seems to be extendable.
4544	!-- Now get the variable ids
4545	i = 1
4546	DO i = 1, dots_num
4547
4548	IF ( statistic_regions == 0 ) THEN
4549	nc_stat = NF90_INQ_VARID( id_set_ts, dots_label(i), &
4550	id_var_dots(i,0) )
4551	CALL netcdf_handle_error( 'netcdf_define_header', 265 )
4552	ELSE
4553	DO j = 0, statistic_regions
4554	WRITE ( suffix, '(''_'',I2.2)' ) j
4555	netcdf_var_name = TRIM( dots_label(i) ) // suffix
4556	nc_stat = NF90_INQ_VARID( id_set_ts, netcdf_var_name, &
4557	id_var_dots(i,j) )
4558	CALL netcdf_handle_error( 'netcdf_define_header', 266 )
4559	ENDDO
4560	ENDIF
4561
4562	ENDDO
4563
4564	!
4565	!-- Update the title attribute on file
4566	!-- In order to avoid 'data mode' errors if updated attributes are larger
4567	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
4568	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
4569	!-- performance loss due to data copying; an alternative strategy would be
4570	!-- to ensure equal attribute size in a job chain. Maybe revise later.
4571	nc_stat = NF90_REDEF( id_set_ts )
4572	CALL netcdf_handle_error( 'netcdf_define_header', 439 )
4573	nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'title', &
4574	TRIM( run_description_header ) )
4575	CALL netcdf_handle_error( 'netcdf_define_header', 267 )
4576	nc_stat = NF90_ENDDEF( id_set_ts )
4577	CALL netcdf_handle_error( 'netcdf_define_header', 440 )
4578	message_string = 'netCDF file for time series ' // &
4579	'from previous run found.' // &
4580	'&This file will be extended.'
4581	CALL message( 'define_netcdf_header', 'PA0259', 0, 0, 0, 6, 0 )
4582
4583
4584	CASE ( 'sp_new' )
4585
4586	!
4587	!-- Define some global attributes of the dataset
4588	IF ( averaging_interval_sp /= 0.0_wp ) THEN
4589	WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
4590	averaging_interval_sp
4591	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
4592	TRIM( run_description_header ) // &
4593	TRIM( time_average_text ) )
4594	CALL netcdf_handle_error( 'netcdf_define_header', 268 )
4595
4596	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval_sp
4597	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'time_avg', &
4598	TRIM( time_average_text ) )
4599	ELSE
4600	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
4601	TRIM( run_description_header ) )
4602	ENDIF
4603	CALL netcdf_handle_error( 'netcdf_define_header', 269 )
4604
4605	!
4606	!-- Define time coordinate for spectra (unlimited dimension)
4607	CALL netcdf_create_dim( id_set_sp, 'time', NF90_UNLIMITED, &
4608	id_dim_time_sp, 270 )
4609	CALL netcdf_create_var( id_set_sp, (/ id_dim_time_sp /), 'time', &
4610	NF90_DOUBLE, id_var_time_sp, 'seconds', '', &
4611	271, 272, 000 )
4612	!
4613	!-- Define the spatial dimensions and coordinates for spectra.
4614	!-- First, determine the number of vertical levels for which spectra
4615	!-- are to be output.
4616	ns = 1
4617	DO WHILE ( comp_spectra_level(ns) /= 999999 .AND. ns <= 100 )
4618	ns = ns + 1
4619	ENDDO
4620	ns = ns - 1
4621
4622	!
4623	!-- Define vertical coordinate grid (zu grid)
4624	CALL netcdf_create_dim( id_set_sp, 'zu_sp', ns, id_dim_zu_sp, 273 )
4625	CALL netcdf_create_var( id_set_sp, (/ id_dim_zu_sp /), 'zu_sp', &
4626	NF90_DOUBLE, id_var_zu_sp, 'meters', '', &
4627	274, 275, 000 )
4628	!
4629	!-- Define vertical coordinate grid (zw grid)
4630	CALL netcdf_create_dim( id_set_sp, 'zw_sp', ns, id_dim_zw_sp, 276 )
4631	CALL netcdf_create_var( id_set_sp, (/ id_dim_zw_sp /), 'zw_sp', &
4632	NF90_DOUBLE, id_var_zw_sp, 'meters', '', &
4633	277, 278, 000 )
4634	!
4635	!-- Define x-axis
4636	CALL netcdf_create_dim( id_set_sp, 'k_x', nx/2, id_dim_x_sp, 279 )
4637	CALL netcdf_create_var( id_set_sp, (/ id_dim_x_sp /), 'k_x', &
4638	NF90_DOUBLE, id_var_x_sp, 'm-1', '', 280, &
4639	281, 000 )
4640	!
4641	!-- Define y-axis
4642	CALL netcdf_create_dim(id_set_sp, 'k_y', ny/2, id_dim_y_sp, 282 )
4643	CALL netcdf_create_var( id_set_sp, (/ id_dim_y_sp /), 'k_y', &
4644	NF90_DOUBLE, id_var_y_sp, 'm-1', '', 283, &
4645	284, 000 )
4646	!
4647	!-- Define the variables
4648	var_list = ';'
4649	i = 1
4650	DO WHILE ( data_output_sp(i) /= ' ' .AND. i <= 10 )
4651	!
4652	!-- First check for the vertical grid
4653	found = .FALSE.
4654	SELECT CASE ( data_output_sp(i) )
4655	!
4656	!-- Most variables are defined on the zu levels
4657	CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', &
4658	'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', &
4659	'rho_ocean', 's', 'sa', 'u', 'v', 'vpt' )
4660
4661	grid_z = 'zu'
4662	!
4663	!-- zw levels
4664	CASE ( 'w' )
4665
4666	grid_z = 'zw'
4667
4668	CASE DEFAULT
4669	!
4670	!-- Check for user-defined quantities (found, grid_x and grid_y
4671	!-- are dummies)
4672	CALL user_define_netcdf_grid( data_output_sp(i), found, &
4673	grid_x, grid_y, grid_z )
4674
4675	END SELECT
4676
4677	IF ( INDEX( spectra_direction(i), 'x' ) /= 0 ) THEN
4678
4679	!
4680	!-- Define the variable
4681	netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
4682	IF ( TRIM( grid_z ) == 'zw' ) THEN
4683	CALL netcdf_create_var( id_set_sp, (/ id_dim_x_sp, &
4684	id_dim_zw_sp, id_dim_time_sp /), &
4685	netcdf_var_name, nc_precision(7), &
4686	id_var_dospx(i), 'unknown', &
4687	netcdf_var_name, 285, 286, 287 )
4688	ELSE
4689	CALL netcdf_create_var( id_set_sp, (/ id_dim_x_sp, &
4690	id_dim_zu_sp, id_dim_time_sp /), &
4691	netcdf_var_name, nc_precision(7), &
4692	id_var_dospx(i), 'unknown', &
4693	netcdf_var_name, 285, 286, 287 )
4694	ENDIF
4695
4696	var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
4697
4698	ENDIF
4699
4700	IF ( INDEX( spectra_direction(i), 'y' ) /= 0 ) THEN
4701
4702	!
4703	!-- Define the variable
4704	netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
4705	IF ( TRIM( grid_z ) == 'zw' ) THEN
4706	CALL netcdf_create_var( id_set_sp, (/ id_dim_y_sp, &
4707	id_dim_zw_sp, id_dim_time_sp /), &
4708	netcdf_var_name, nc_precision(7), &
4709	id_var_dospy(i), 'unknown', &
4710	netcdf_var_name, 288, 289, 290 )
4711	ELSE
4712	CALL netcdf_create_var( id_set_sp, (/ id_dim_y_sp, &
4713	id_dim_zu_sp, id_dim_time_sp /), &
4714	netcdf_var_name, nc_precision(7), &
4715	id_var_dospy(i), 'unknown', &
4716	netcdf_var_name, 288, 289, 290 )
4717	ENDIF
4718
4719	var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
4720
4721	ENDIF
4722
4723	i = i + 1
4724
4725	ENDDO
4726
4727	!
4728	!-- Write the list of variables as global attribute (this is used by
4729	!-- restart runs)
4730	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', var_list )
4731	CALL netcdf_handle_error( 'netcdf_define_header', 291 )
4732
4733	!
4734	!-- Leave netCDF define mode
4735	nc_stat = NF90_ENDDEF( id_set_sp )
4736	CALL netcdf_handle_error( 'netcdf_define_header', 292 )
4737
4738	!
4739	!-- Write axis data: zu_sp, zw_sp, k_x, k_y
4740	ALLOCATE( netcdf_data(1:ns) )
4741
4742	!
4743	!-- Write zu data
4744	netcdf_data(1:ns) = zu( comp_spectra_level(1:ns) )
4745	nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zu_sp, netcdf_data, &
4746	start = (/ 1 /), count = (/ ns /) )
4747	CALL netcdf_handle_error( 'netcdf_define_header', 293 )
4748
4749	!
4750	!-- Write zw data
4751	netcdf_data(1:ns) = zw( comp_spectra_level(1:ns) )
4752	nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zw_sp, netcdf_data, &
4753	start = (/ 1 /), count = (/ ns /) )
4754	CALL netcdf_handle_error( 'netcdf_define_header', 294 )
4755
4756	DEALLOCATE( netcdf_data )
4757
4758	!
4759	!-- Write data for x and y axis (wavenumbers)
4760	ALLOCATE( netcdf_data(nx/2) )
4761	DO i = 1, nx/2
4762	netcdf_data(i) = 2.0_wp * pi * i / ( dx * ( nx + 1 ) )
4763	ENDDO
4764
4765	nc_stat = NF90_PUT_VAR( id_set_sp, id_var_x_sp, netcdf_data, &
4766	start = (/ 1 /), count = (/ nx/2 /) )
4767	CALL netcdf_handle_error( 'netcdf_define_header', 295 )
4768
4769	DEALLOCATE( netcdf_data )
4770
4771	ALLOCATE( netcdf_data(ny/2) )
4772	DO i = 1, ny/2
4773	netcdf_data(i) = 2.0_wp * pi * i / ( dy * ( ny + 1 ) )
4774	ENDDO
4775
4776	nc_stat = NF90_PUT_VAR( id_set_sp, id_var_y_sp, netcdf_data, &
4777	start = (/ 1 /), count = (/ ny/2 /) )
4778	CALL netcdf_handle_error( 'netcdf_define_header', 296 )
4779
4780	DEALLOCATE( netcdf_data )
4781
4782
4783	CASE ( 'sp_ext' )
4784
4785	!
4786	!-- Get the list of variables and compare with the actual run.
4787	!-- First var_list_old has to be reset, since GET_ATT does not assign
4788	!-- trailing blanks.
4789	var_list_old = ' '
4790	nc_stat = NF90_GET_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', &
4791	var_list_old )
4792	CALL netcdf_handle_error( 'netcdf_define_header', 297 )
4793	var_list = ';'
4794	i = 1
4795	DO WHILE ( data_output_sp(i) /= ' ' .AND. i <= 10 )
4796
4797	IF ( INDEX( spectra_direction(i), 'x' ) /= 0 ) THEN
4798	netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
4799	var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
4800	ENDIF
4801
4802	IF ( INDEX( spectra_direction(i), 'y' ) /= 0 ) THEN
4803	netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
4804	var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
4805	ENDIF
4806
4807	i = i + 1
4808
4809	ENDDO
4810
4811	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
4812	message_string = 'netCDF file for spectra ' // &
4813	'from previous run found,' // &
4814	'& but this file cannot be extended due to' // &
4815	' variable mismatch.' // &
4816	'&New file is created instead.'
4817	CALL message( 'define_netcdf_header', 'PA0260', 0, 1, 0, 6, 0 )
4818	extend = .FALSE.
4819	RETURN
4820	ENDIF
4821
4822	!
4823	!-- Determine the number of current vertical levels for which spectra
4824	!-- shall be output
4825	ns = 1
4826	DO WHILE ( comp_spectra_level(ns) /= 999999 .AND. ns <= 100 )
4827	ns = ns + 1
4828	ENDDO
4829	ns = ns - 1
4830
4831	!
4832	!-- Get and compare the number of vertical levels
4833	nc_stat = NF90_INQ_VARID( id_set_sp, 'zu_sp', id_var_zu_sp )
4834	CALL netcdf_handle_error( 'netcdf_define_header', 298 )
4835
4836	nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_zu_sp, &
4837	dimids = id_dim_zu_sp_old )
4838	CALL netcdf_handle_error( 'netcdf_define_header', 299 )
4839	id_dim_zu_sp = id_dim_zu_sp_old(1)
4840
4841	nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_zu_sp, &
4842	len = ns_old )
4843	CALL netcdf_handle_error( 'netcdf_define_header', 300 )
4844
4845	IF ( ns /= ns_old ) THEN
4846	message_string = 'netCDF file for spectra ' // &
4847	' from previous run found,' // &
4848	'&but this file cannot be extended due to' // &
4849	' mismatch in number of' // &
4850	'&vertical levels.' // &
4851	'&New file is created instead.'
4852	CALL message( 'define_netcdf_header', 'PA0261', 0, 1, 0, 6, 0 )
4853	extend = .FALSE.
4854	RETURN
4855	ENDIF
4856
4857	!
4858	!-- Get and compare the heights of the cross sections
4859	ALLOCATE( netcdf_data(1:ns_old) )
4860
4861	nc_stat = NF90_GET_VAR( id_set_sp, id_var_zu_sp, netcdf_data )
4862	CALL netcdf_handle_error( 'netcdf_define_header', 301 )
4863
4864	DO i = 1, ns
4865	IF ( zu(comp_spectra_level(i)) /= netcdf_data(i) ) THEN
4866	message_string = 'netCDF file for spectra ' // &
4867	' from previous run found,' // &
4868	'&but this file cannot be extended due to' // &
4869	' mismatch in heights of' // &
4870	'&vertical levels.' // &
4871	'&New file is created instead.'
4872	CALL message( 'define_netcdf_header', 'PA0262', 0, 1, 0, 6, 0 )
4873	extend = .FALSE.
4874	RETURN
4875	ENDIF
4876	ENDDO
4877
4878	DEALLOCATE( netcdf_data )
4879
4880	!
4881	!-- Get the id of the time coordinate (unlimited coordinate) and its
4882	!-- last index on the file. The next time level is plsp..count+1.
4883	!-- The current time must be larger than the last output time
4884	!-- on the file.
4885	nc_stat = NF90_INQ_VARID( id_set_sp, 'time', id_var_time_sp )
4886	CALL netcdf_handle_error( 'netcdf_define_header', 302 )
4887
4888	nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_time_sp, &
4889	dimids = id_dim_time_old )
4890	CALL netcdf_handle_error( 'netcdf_define_header', 303 )
4891	id_dim_time_sp = id_dim_time_old(1)
4892
4893	nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_time_sp, &
4894	len = dosp_time_count )
4895	CALL netcdf_handle_error( 'netcdf_define_header', 304 )
4896
4897	nc_stat = NF90_GET_VAR( id_set_sp, id_var_time_sp, &
4898	last_time_coordinate, &
4899	start = (/ dosp_time_count /), &
4900	count = (/ 1 /) )
4901	CALL netcdf_handle_error( 'netcdf_define_header', 305 )
4902
4903	IF ( last_time_coordinate(1) >= simulated_time ) THEN
4904	message_string = 'netCDF file for spectra ' // &
4905	'from previous run found,' // &
4906	'&but this file cannot be extended becaus' // &
4907	'e the current output time' // &
4908	'&is less or equal than the last output t' // &
4909	'ime on this file.' // &
4910	'&New file is created instead.'
4911	CALL message( 'define_netcdf_header', 'PA0263', 0, 1, 0, 6, 0 )
4912	dosp_time_count = 0
4913	extend = .FALSE.
4914	RETURN
4915	ENDIF
4916
4917	!
4918	!-- Dataset seems to be extendable.
4919	!-- Now get the variable ids.
4920	i = 1
4921	DO WHILE ( data_output_sp(i) /= ' ' .AND. i <= 10 )
4922
4923	IF ( INDEX( spectra_direction(i), 'x' ) /= 0 ) THEN
4924	netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
4925	nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
4926	id_var_dospx(i) )
4927	CALL netcdf_handle_error( 'netcdf_define_header', 306 )
4928	ENDIF
4929
4930	IF ( INDEX( spectra_direction(i), 'y' ) /= 0 ) THEN
4931	netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
4932	nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
4933	id_var_dospy(i) )
4934	CALL netcdf_handle_error( 'netcdf_define_header', 307 )
4935	ENDIF
4936
4937	i = i + 1
4938
4939	ENDDO
4940
4941	!
4942	!-- Update the title attribute on file
4943	!-- In order to avoid 'data mode' errors if updated attributes are larger
4944	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
4945	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
4946	!-- performance loss due to data copying; an alternative strategy would be
4947	!-- to ensure equal attribute size in a job chain. Maybe revise later.
4948	nc_stat = NF90_REDEF( id_set_sp )
4949	CALL netcdf_handle_error( 'netcdf_define_header', 441 )
4950	IF ( averaging_interval_sp /= 0.0_wp ) THEN
4951	WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
4952	averaging_interval_sp
4953	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
4954	TRIM( run_description_header ) // &
4955	TRIM( time_average_text ) )
4956	CALL netcdf_handle_error( 'netcdf_define_header', 308 )
4957
4958	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval_sp
4959	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'time_avg', &
4960	TRIM( time_average_text ) )
4961	ELSE
4962	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
4963	TRIM( run_description_header ) )
4964	ENDIF
4965	CALL netcdf_handle_error( 'netcdf_define_header', 309 )
4966	nc_stat = NF90_ENDDEF( id_set_sp )
4967	CALL netcdf_handle_error( 'netcdf_define_header', 442 )
4968	message_string = 'netCDF file for spectra ' // &
4969	'from previous run found.' // &
4970	'&This file will be extended.'
4971	CALL message( 'define_netcdf_header', 'PA0264', 0, 0, 0, 6, 0 )
4972
4973
4974	CASE ( 'pt_new' )
4975
4976	!
4977	!-- Define some global attributes of the dataset
4978	nc_stat = NF90_PUT_ATT( id_set_prt, NF90_GLOBAL, 'title', &
4979	TRIM( run_description_header ) )
4980	CALL netcdf_handle_error( 'netcdf_define_header', 310 )
4981
4982	!
4983	!-- Define time coordinate for particles (unlimited dimension)
4984	CALL netcdf_create_dim( id_set_prt, 'time', NF90_UNLIMITED, &
4985	id_dim_time_prt, 311 )
4986	CALL netcdf_create_var( id_set_prt, (/ id_dim_time_prt /), 'time', &
4987	NF90_DOUBLE, id_var_time_prt, 'seconds', '', &
4988	312, 313, 000 )
4989	!
4990	!-- netCDF4 allows more than one unlimited dimension
4991	CALL netcdf_create_dim( id_set_prt, 'particle_number', &
4992	NF90_UNLIMITED, id_dim_prtnum, 314 )
4993
4994	CALL netcdf_create_var( id_set_prt, (/ id_dim_prtnum /), &
4995	'particle_number', NF90_DOUBLE, &
4996	id_var_prtnum, 'particle number', '', 315, &
4997	316, 000 )
4998	!
4999	!-- Define variable which contains the real number of particles in use
5000	CALL netcdf_create_var( id_set_prt, (/ id_dim_time_prt /), &
5001	'real_num_of_prt', NF90_DOUBLE, &
5002	id_var_rnop_prt, 'particle number', '', 317, &
5003	318, 000 )
5004	!
5005	!-- Define the variables
5006	DO i = 1, 17
5007
5008	CALL netcdf_create_var( id_set_prt, (/ id_dim_prtnum, &
5009	id_dim_time_prt /), prt_var_names(i), &
5010	nc_precision(8), id_var_prt(i), &
5011	TRIM( prt_var_units(i) ), &
5012	TRIM( prt_var_names(i) ), 319, 320, 321 )
5013
5014	ENDDO
5015
5016	!
5017	!-- Leave netCDF define mode
5018	nc_stat = NF90_ENDDEF( id_set_prt )
5019	CALL netcdf_handle_error( 'netcdf_define_header', 322 )
5020
5021
5022	CASE ( 'pt_ext' )
5023
5024	!
5025	!-- Get the id of the time coordinate (unlimited coordinate) and its
5026	!-- last index on the file. The next time level is prt..count+1.
5027	!-- The current time must be larger than the last output time
5028	!-- on the file.
5029	nc_stat = NF90_INQ_VARID( id_set_prt, 'time', id_var_time_prt )
5030	CALL netcdf_handle_error( 'netcdf_define_header', 323 )
5031
5032	nc_stat = NF90_INQUIRE_VARIABLE( id_set_prt, id_var_time_prt, &
5033	dimids = id_dim_time_old )
5034	CALL netcdf_handle_error( 'netcdf_define_header', 324 )
5035	id_dim_time_prt = id_dim_time_old(1)
5036
5037	nc_stat = NF90_INQUIRE_DIMENSION( id_set_prt, id_dim_time_prt, &
5038	len = prt_time_count )
5039	CALL netcdf_handle_error( 'netcdf_define_header', 325 )
5040
5041	nc_stat = NF90_GET_VAR( id_set_prt, id_var_time_prt, &
5042	last_time_coordinate, &
5043	start = (/ prt_time_count /), &
5044	count = (/ 1 /) )
5045	CALL netcdf_handle_error( 'netcdf_define_header', 326 )
5046
5047	IF ( last_time_coordinate(1) >= simulated_time ) THEN
5048	message_string = 'netCDF file for particles ' // &
5049	'from previous run found,' // &
5050	'&but this file cannot be extended becaus' // &
5051	'e the current output time' // &
5052	'&is less or equal than the last output t' // &
5053	'ime on this file.' // &
5054	'&New file is created instead.'
5055	CALL message( 'define_netcdf_header', 'PA0265', 0, 1, 0, 6, 0 )
5056	prt_time_count = 0
5057	extend = .FALSE.
5058	RETURN
5059	ENDIF
5060
5061	!
5062	!-- Dataset seems to be extendable.
5063	!-- Now get the variable ids.
5064	nc_stat = NF90_INQ_VARID( id_set_prt, 'real_num_of_prt', &
5065	id_var_rnop_prt )
5066	CALL netcdf_handle_error( 'netcdf_define_header', 327 )
5067
5068	DO i = 1, 17
5069
5070	nc_stat = NF90_INQ_VARID( id_set_prt, prt_var_names(i), &
5071	id_var_prt(i) )
5072	CALL netcdf_handle_error( 'netcdf_define_header', 328 )
5073
5074	ENDDO
5075
5076	message_string = 'netCDF file for particles ' // &
5077	'from previous run found.' // &
5078	'&This file will be extended.'
5079	CALL message( 'define_netcdf_header', 'PA0266', 0, 0, 0, 6, 0 )
5080
5081
5082
5083	CASE ( 'ps_new' )
5084
5085	!
5086	!-- Define some global attributes of the dataset
5087	nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
5088	TRIM( run_description_header ) )
5089	CALL netcdf_handle_error( 'netcdf_define_header', 396 )
5090
5091	!
5092	!-- Define time coordinate for particle time series (unlimited dimension)
5093	CALL netcdf_create_dim( id_set_pts, 'time', NF90_UNLIMITED, &
5094	id_dim_time_pts, 397 )
5095	CALL netcdf_create_var( id_set_pts, (/ id_dim_time_pts /), 'time', &
5096	NF90_DOUBLE, id_var_time_pts, 'seconds', '', &
5097	398, 399, 000 )
5098	!
5099	!-- Define the variables. If more than one particle group is defined,
5100	!-- define seperate variables for each group
5101	var_list = ';'
5102	DO i = 1, dopts_num
5103
5104	DO j = 0, number_of_particle_groups
5105
5106	IF ( j == 0 ) THEN
5107	suffix = ''
5108	ELSE
5109	WRITE ( suffix, '(''_'',I2.2)' ) j
5110	ENDIF
5111
5112	IF ( j == 0 ) THEN
5113	CALL netcdf_create_var( id_set_pts, (/ id_dim_time_pts /), &
5114	TRIM( dopts_label(i) ) // suffix, &
5115	nc_precision(6), id_var_dopts(i,j), &
5116	TRIM( dopts_unit(i) ), &
5117	TRIM( dopts_label(i) ), 400, 401, &
5118	402 )
5119	ELSE
5120	CALL netcdf_create_var( id_set_pts, (/ id_dim_time_pts /), &
5121	TRIM( dopts_label(i) ) // suffix, &
5122	nc_precision(6), id_var_dopts(i,j), &
5123	TRIM( dopts_unit(i) ), &
5124	TRIM( dopts_label(i) ) // ' PG ' // &
5125	suffix(2:3), 400, 401, 402 )
5126	ENDIF
5127
5128	var_list = TRIM( var_list ) // TRIM( dopts_label(i) ) // &
5129	suffix // ';'
5130
5131	IF ( number_of_particle_groups == 1 ) EXIT
5132
5133	ENDDO
5134
5135	ENDDO
5136
5137	!
5138	!-- Write the list of variables as global attribute (this is used by
5139	!-- restart runs)
5140	nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
5141	var_list )
5142	CALL netcdf_handle_error( 'netcdf_define_header', 403 )
5143
5144
5145	!
5146	!-- Leave netCDF define mode
5147	nc_stat = NF90_ENDDEF( id_set_pts )
5148	CALL netcdf_handle_error( 'netcdf_define_header', 404 )
5149
5150
5151	CASE ( 'ps_ext' )
5152
5153	!
5154	!-- Get the list of variables and compare with the actual run.
5155	!-- First var_list_old has to be reset, since GET_ATT does not assign
5156	!-- trailing blanks.
5157	var_list_old = ' '
5158	nc_stat = NF90_GET_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
5159	var_list_old )
5160	CALL netcdf_handle_error( 'netcdf_define_header', 405 )
5161
5162	var_list = ';'
5163	i = 1
5164	DO i = 1, dopts_num
5165
5166	DO j = 0, number_of_particle_groups
5167
5168	IF ( j == 0 ) THEN
5169	suffix = ''
5170	ELSE
5171	WRITE ( suffix, '(''_'',I2.2)' ) j
5172	ENDIF
5173
5174	var_list = TRIM( var_list ) // TRIM( dopts_label(i) ) // &
5175	suffix // ';'
5176
5177	IF ( number_of_particle_groups == 1 ) EXIT
5178
5179	ENDDO
5180
5181	ENDDO
5182
5183	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
5184	message_string = 'netCDF file for particle time series ' // &
5185	'from previous run found,' // &
5186	'& but this file cannot be extended due to' // &
5187	' variable mismatch.' // &
5188	'&New file is created instead.'
5189	CALL message( 'define_netcdf_header', 'PA0267', 0, 1, 0, 6, 0 )
5190	extend = .FALSE.
5191	RETURN
5192	ENDIF
5193
5194	!
5195	!-- Get the id of the time coordinate (unlimited coordinate) and its
5196	!-- last index on the file. The next time level is dots..count+1.
5197	!-- The current time must be larger than the last output time
5198	!-- on the file.
5199	nc_stat = NF90_INQ_VARID( id_set_pts, 'time', id_var_time_pts )
5200	CALL netcdf_handle_error( 'netcdf_define_header', 406 )
5201
5202	nc_stat = NF90_INQUIRE_VARIABLE( id_set_pts, id_var_time_pts, &
5203	dimids = id_dim_time_old )
5204	CALL netcdf_handle_error( 'netcdf_define_header', 407 )
5205	id_dim_time_pts = id_dim_time_old(1)
5206
5207	nc_stat = NF90_INQUIRE_DIMENSION( id_set_pts, id_dim_time_pts, &
5208	len = dopts_time_count )
5209	CALL netcdf_handle_error( 'netcdf_define_header', 408 )
5210
5211	nc_stat = NF90_GET_VAR( id_set_pts, id_var_time_pts, &
5212	last_time_coordinate, &
5213	start = (/ dopts_time_count /), &
5214	count = (/ 1 /) )
5215	CALL netcdf_handle_error( 'netcdf_define_header', 409 )
5216
5217	IF ( last_time_coordinate(1) >= simulated_time ) THEN
5218	message_string = 'netCDF file for particle time series ' // &
5219	'from previous run found,' // &
5220	'&but this file cannot be extended becaus' // &
5221	'e the current output time' // &
5222	'&is less or equal than the last output t' // &
5223	'ime on this file.' // &
5224	'&New file is created instead.'
5225	CALL message( 'define_netcdf_header', 'PA0268', 0, 1, 0, 6, 0 )
5226	dopts_time_count = 0
5227	extend = .FALSE.
5228	RETURN
5229	ENDIF
5230
5231	!
5232	!-- Dataset seems to be extendable.
5233	!-- Now get the variable ids
5234	i = 1
5235	DO i = 1, dopts_num
5236
5237	DO j = 0, number_of_particle_groups
5238
5239	IF ( j == 0 ) THEN
5240	suffix = ''
5241	ELSE
5242	WRITE ( suffix, '(''_'',I2.2)' ) j
5243	ENDIF
5244
5245	netcdf_var_name = TRIM( dopts_label(i) ) // suffix
5246
5247	nc_stat = NF90_INQ_VARID( id_set_pts, netcdf_var_name, &
5248	id_var_dopts(i,j) )
5249	CALL netcdf_handle_error( 'netcdf_define_header', 410 )
5250
5251	IF ( number_of_particle_groups == 1 ) EXIT
5252
5253	ENDDO
5254
5255	ENDDO
5256
5257	!
5258	!-- Update the title attribute on file
5259	!-- In order to avoid 'data mode' errors if updated attributes are larger
5260	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
5261	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
5262	!-- performance loss due to data copying; an alternative strategy would be
5263	!-- to ensure equal attribute size in a job chain. Maybe revise later.
5264	nc_stat = NF90_REDEF( id_set_pts )
5265	CALL netcdf_handle_error( 'netcdf_define_header', 443 )
5266	nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
5267	TRIM( run_description_header ) )
5268	CALL netcdf_handle_error( 'netcdf_define_header', 411 )
5269	nc_stat = NF90_ENDDEF( id_set_pts )
5270	CALL netcdf_handle_error( 'netcdf_define_header', 444 )
5271	message_string = 'netCDF file for particle time series ' // &
5272	'from previous run found.' // &
5273	'&This file will be extended.'
5274	CALL message( 'netcdf_define_header', 'PA0269', 0, 0, 0, 6, 0 )
5275
5276	!
5277	!-- Flight data
5278	CASE ( 'fl_new' )
5279	!
5280	!-- Define some global attributes of the dataset
5281	nc_stat = NF90_PUT_ATT( id_set_fl, NF90_GLOBAL, 'title', &
5282	TRIM( run_description_header ) )
5283	CALL netcdf_handle_error( 'netcdf_define_header', 249 )
5284
5285	!
5286	!-- Define time and location coordinates for flight space-time series
5287	!-- (unlimited dimension)
5288	!-- Error number must still be set appropriately.
5289	CALL netcdf_create_dim( id_set_fl, 'time', NF90_UNLIMITED, &
5290	id_dim_time_fl, 250 )
5291	CALL netcdf_create_var( id_set_fl, (/ id_dim_time_fl /), 'time', &
5292	NF90_DOUBLE, id_var_time_fl, 'seconds', '', &
5293	251, 252, 000 )
5294
5295	DO l = 1, num_leg
5296	CALL netcdf_create_dim( id_set_fl, dofl_dim_label_x(l), &
5297	NF90_UNLIMITED, id_dim_x_fl(l), 250 )
5298	CALL netcdf_create_dim( id_set_fl, dofl_dim_label_y(l), &
5299	NF90_UNLIMITED, id_dim_y_fl(l), 250 )
5300	CALL netcdf_create_dim( id_set_fl, dofl_dim_label_z(l), &
5301	NF90_UNLIMITED, id_dim_z_fl(l), 250 )
5302
5303	CALL netcdf_create_var( id_set_fl, (/ id_dim_x_fl(l) /), &
5304	dofl_dim_label_x(l), NF90_DOUBLE, &
5305	id_var_x_fl(l), 'm', '', 251, 252, 000 )
5306	CALL netcdf_create_var( id_set_fl, (/ id_dim_y_fl(l) /), &
5307	dofl_dim_label_y(l), NF90_DOUBLE, &
5308	id_var_y_fl(l), 'm', '', 251, 252, 000 )
5309	CALL netcdf_create_var( id_set_fl, (/ id_dim_z_fl(l) /), &
5310	dofl_dim_label_z(l), NF90_DOUBLE, &
5311	id_var_z_fl(l), 'm', '', 251, 252, 000 )
5312	ENDDO
5313	!
5314	!-- Define the variables
5315	var_list = ';'
5316	k = 1
5317	DO l = 1, num_leg
5318	DO i = 1, num_var_fl
5319
5320	CALL netcdf_create_var( id_set_fl, (/ id_dim_time_fl /), &
5321	dofl_label(k), nc_precision(9), &
5322	id_var_dofl(k), &
5323	TRIM( dofl_unit(k) ), &
5324	TRIM( dofl_label(k) ), 253, 254, 255 )
5325
5326	k = k + 1
5327
5328	ENDDO
5329
5330	ENDDO
5331
5332	!
5333	!-- Write the list of variables as global attribute (this is used by
5334	!-- restart runs)
5335	nc_stat = NF90_PUT_ATT( id_set_fl, NF90_GLOBAL, 'VAR_LIST', var_list )
5336	CALL netcdf_handle_error( 'netcdf_define_header', 258 )
5337
5338	!
5339	!-- Leave netCDF define mode
5340	nc_stat = NF90_ENDDEF( id_set_fl )
5341	CALL netcdf_handle_error( 'netcdf_define_header', 259 )
5342
5343
5344	CASE ( 'fl_ext' )
5345
5346	!
5347	!-- Get the list of variables and compare with the actual run.
5348	!-- First var_list_old has to be reset, since GET_ATT does not assign
5349	!-- trailing blanks.
5350	var_list_old = ' '
5351	nc_stat = NF90_GET_ATT( id_set_fl, NF90_GLOBAL, 'VAR_LIST', &
5352	var_list_old )
5353	CALL netcdf_handle_error( 'netcdf_define_header', 260 )
5354
5355	var_list = ';'
5356	i = 1
5357	DO i = 1, num_leg * num_var_fl
5358
5359	var_list = TRIM( var_list ) // TRIM( dofl_label(i) ) // ';'
5360
5361	ENDDO
5362
5363	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
5364	message_string = 'netCDF file for flight time series ' // &
5365	'from previous run found,' // &
5366	'& but this file cannot be extended due to' // &
5367	' variable mismatch.' // &
5368	'&New file is created instead.'
5369	CALL message( 'define_netcdf_header', 'PA0257', 0, 1, 0, 6, 0 )
5370	extend = .FALSE.
5371	RETURN
5372	ENDIF
5373
5374	!
5375	!-- Get the id of the time coordinate (unlimited coordinate) and its
5376	!-- last index on the file. The next time level is dofl_time_count+1.
5377	!-- The current time must be larger than the last output time
5378	!-- on the file.
5379	nc_stat = NF90_INQ_VARID( id_set_fl, 'time', id_var_time_fl )
5380	CALL netcdf_handle_error( 'netcdf_define_header', 261 )
5381
5382	nc_stat = NF90_INQUIRE_VARIABLE( id_set_fl, id_var_time_fl, &
5383	dimids = id_dim_time_old )
5384	CALL netcdf_handle_error( 'netcdf_define_header', 262 )
5385	id_dim_time_fl = id_dim_time_old(1)
5386
5387	nc_stat = NF90_INQUIRE_DIMENSION( id_set_fl, id_dim_time_fl, &
5388	len = dofl_time_count )
5389	CALL netcdf_handle_error( 'netcdf_define_header', 263 )
5390
5391	nc_stat = NF90_GET_VAR( id_set_fl, id_var_time_fl, &
5392	last_time_coordinate, &
5393	start = (/ dofl_time_count /), &
5394	count = (/ 1 /) )
5395	CALL netcdf_handle_error( 'netcdf_define_header', 264 )
5396
5397	IF ( last_time_coordinate(1) >= simulated_time ) THEN
5398	message_string = 'netCDF file for flight-time series ' // &
5399	'from previous run found,' // &
5400	'&but this file cannot be extended becaus' // &
5401	'e the current output time' // &
5402	'&is less or equal than the last output t' // &
5403	'ime on this file.' // &
5404	'&New file is created instead.'
5405	CALL message( 'define_netcdf_header', 'PA0258', 0, 1, 0, 6, 0 )
5406	dofl_time_count = 0
5407	extend = .FALSE.
5408	RETURN
5409	ENDIF
5410
5411	!
5412	!-- Dataset seems to be extendable.
5413	!-- Now get the remaining dimension and variable ids
5414	DO l = 1, num_leg
5415	nc_stat = NF90_INQ_VARID( id_set_fl, dofl_dim_label_x(l), &
5416	id_var_x_fl(l) )
5417	CALL netcdf_handle_error( 'netcdf_define_header', 265 )
5418	nc_stat = NF90_INQ_VARID( id_set_fl, dofl_dim_label_y(l), &
5419	id_var_y_fl(l) )
5420	CALL netcdf_handle_error( 'netcdf_define_header', 265 )
5421	nc_stat = NF90_INQ_VARID( id_set_fl, dofl_dim_label_z(l), &
5422	id_var_z_fl(l) )
5423	CALL netcdf_handle_error( 'netcdf_define_header', 265 )
5424
5425	ENDDO
5426
5427
5428	DO i = 1, num_leg * num_var_fl
5429
5430	nc_stat = NF90_INQ_VARID( id_set_fl, dofl_label(i), &
5431	id_var_dofl(i) )
5432	CALL netcdf_handle_error( 'netcdf_define_header', 265 )
5433
5434	ENDDO
5435
5436	!
5437	!-- Update the title attribute on file
5438	!-- In order to avoid 'data mode' errors if updated attributes are larger
5439	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
5440	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
5441	!-- performance loss due to data copying; an alternative strategy would be
5442	!-- to ensure equal attribute size in a job chain. Maybe revise later.
5443	nc_stat = NF90_REDEF( id_set_fl )
5444	CALL netcdf_handle_error( 'netcdf_define_header', 439 )
5445	nc_stat = NF90_PUT_ATT( id_set_fl, NF90_GLOBAL, 'title', &
5446	TRIM( run_description_header ) )
5447	CALL netcdf_handle_error( 'netcdf_define_header', 267 )
5448	nc_stat = NF90_ENDDEF( id_set_fl )
5449	CALL netcdf_handle_error( 'netcdf_define_header', 440 )
5450	message_string = 'netCDF file for flight-time series ' // &
5451	'from previous run found.' // &
5452	'&This file will be extended.'
5453	CALL message( 'define_netcdf_header', 'PA0259', 0, 0, 0, 6, 0 )
5454
5455
5456	CASE DEFAULT
5457
5458	message_string = 'mode "' // TRIM( mode) // '" not supported'
5459	CALL message( 'netcdf_define_header', 'PA0270', 0, 0, 0, 6, 0 )
5460
5461	END SELECT
5462
5463	#endif
5464	END SUBROUTINE netcdf_define_header
5465
5466
5467	!------------------------------------------------------------------------------!
5468	! Description:
5469	! ------------
5470	!> Creates a netCDF file and give back the id. The parallel flag has to be TRUE
5471	!> for parallel netCDF output support.
5472	!------------------------------------------------------------------------------!
5473
5474	SUBROUTINE netcdf_create_file( filename , id, parallel, errno )
5475	#if defined( __netcdf )
5476
5477	USE pegrid
5478
5479	IMPLICIT NONE
5480
5481	CHARACTER (LEN=*), INTENT(IN) :: filename
5482	INTEGER, INTENT(IN) :: errno
5483	INTEGER, INTENT(OUT) :: id
5484	LOGICAL, INTENT(IN) :: parallel
5485
5486
5487	!
5488	!-- Create a new netCDF output file with requested netCDF format
5489	IF ( netcdf_data_format == 1 ) THEN
5490	!
5491	!-- Classic netCDF format
5492	nc_stat = NF90_CREATE( filename, NF90_NOCLOBBER, id )
5493
5494	ELSEIF ( netcdf_data_format == 2 ) THEN
5495	!
5496	!-- 64bit-offset format
5497	nc_stat = NF90_CREATE( filename, &
5498	IOR( NF90_NOCLOBBER, NF90_64BIT_OFFSET ), id )
5499
5500	#if defined( __netcdf4 )
5501	ELSEIF ( netcdf_data_format == 3 .OR. &
5502	( .NOT. parallel .AND. netcdf_data_format == 5 ) ) THEN
5503	!
5504	!-- netCDF4/HDF5 format
5505	nc_stat = NF90_CREATE( filename, IOR( NF90_NOCLOBBER, NF90_NETCDF4 ), id )
5506
5507	ELSEIF ( netcdf_data_format == 4 .OR. &
5508	( .NOT. parallel .AND. netcdf_data_format == 6 ) ) THEN
5509	!
5510	!-- netCDF4/HDF5 format with classic model flag
5511	nc_stat = NF90_CREATE( filename, &
5512	IOR( NF90_NOCLOBBER, &
5513	IOR( NF90_CLASSIC_MODEL, NF90_HDF5 ) ), id )
5514
5515	#if defined( __netcdf4_parallel )
5516	ELSEIF ( netcdf_data_format == 5 .AND. parallel ) THEN
5517	!
5518	!-- netCDF4/HDF5 format, parallel
5519	nc_stat = NF90_CREATE( filename, &
5520	IOR( NF90_NOCLOBBER, &
5521	IOR( NF90_NETCDF4, NF90_MPIIO ) ), &
5522	id, COMM = comm2d, INFO = MPI_INFO_NULL )
5523
5524	ELSEIF ( netcdf_data_format == 6 .AND. parallel ) THEN
5525	!
5526	!-- netCDF4/HDF5 format with classic model flag, parallel
5527	nc_stat = NF90_CREATE( filename, &
5528	IOR( NF90_NOCLOBBER, &
5529	IOR( NF90_MPIIO, &
5530	IOR( NF90_CLASSIC_MODEL, NF90_HDF5 ) ) ), &
5531	id, COMM = comm2d, INFO = MPI_INFO_NULL )
5532
5533	#endif
5534	#endif
5535	ENDIF
5536
5537	CALL netcdf_handle_error( 'netcdf_create_file', errno )
5538	#endif
5539	END SUBROUTINE netcdf_create_file
5540
5541	!------------------------------------------------------------------------------!
5542	! Description:
5543	! ------------
5544	!> Opens an existing netCDF file for writing and gives back the id.
5545	!> The parallel flag has to be TRUE for parallel netCDF output support.
5546	!------------------------------------------------------------------------------!
5547	SUBROUTINE netcdf_open_write_file( filename, id, parallel, errno )
5548	#if defined( __netcdf )
5549
5550	USE pegrid
5551
5552	IMPLICIT NONE
5553
5554	CHARACTER (LEN=*), INTENT(IN) :: filename
5555	INTEGER, INTENT(IN) :: errno
5556	INTEGER, INTENT(OUT) :: id
5557	LOGICAL, INTENT(IN) :: parallel
5558
5559
5560	IF ( netcdf_data_format < 5 .OR. .NOT. parallel ) THEN
5561	nc_stat = NF90_OPEN( filename, NF90_WRITE, id )
5562	#if defined( __netcdf4 )
5563	#if defined( __netcdf4_parallel )
5564	ELSEIF ( netcdf_data_format > 4 .AND. parallel ) THEN
5565	nc_stat = NF90_OPEN( filename, IOR( NF90_WRITE, NF90_MPIIO ), id, &
5566	COMM = comm2d, INFO = MPI_INFO_NULL )
5567	#endif
5568	#endif
5569	ENDIF
5570
5571	CALL netcdf_handle_error( 'netcdf_open_write_file', errno )
5572	#endif
5573	END SUBROUTINE netcdf_open_write_file
5574
5575
5576	!------------------------------------------------------------------------------!
5577	! Description:
5578	! ------------
5579	!> Prints out a text message corresponding to the current status.
5580	!------------------------------------------------------------------------------!
5581
5582	SUBROUTINE netcdf_handle_error( routine_name, errno )
5583	#if defined( __netcdf )
5584
5585
5586	USE control_parameters, &
5587	ONLY: message_string
5588
5589	IMPLICIT NONE
5590
5591	CHARACTER(LEN=6) :: message_identifier
5592	CHARACTER(LEN=*) :: routine_name
5593
5594	INTEGER(iwp) :: errno
5595
5596	IF ( nc_stat /= NF90_NOERR ) THEN
5597
5598	WRITE( message_identifier, '(''NC'',I4.4)' ) errno
5599	message_string = TRIM( NF90_STRERROR( nc_stat ) )
5600
5601	CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
5602
5603	ENDIF
5604
5605	#endif
5606	END SUBROUTINE netcdf_handle_error
5607
5608
5609	!------------------------------------------------------------------------------!
5610	! Description:
5611	! ------------
5612	!> Create a dimension in NetCDF file
5613	!------------------------------------------------------------------------------!
5614
5615	SUBROUTINE netcdf_create_dim(ncid, dim_name, ncdim_type, ncdim_id, error_no)
5616
5617	#if defined( __netcdf )
5618
5619	USE kinds
5620
5621	IMPLICIT NONE
5622
5623	CHARACTER(LEN=*), INTENT(IN) :: dim_name
5624
5625	INTEGER, INTENT(IN) :: error_no
5626	INTEGER, INTENT(IN) :: ncid
5627	INTEGER, INTENT(OUT) :: ncdim_id
5628	INTEGER, INTENT(IN) :: ncdim_type
5629
5630	!
5631	!-- Define time coordinate for volume data (unlimited dimension)
5632	nc_stat = NF90_DEF_DIM( ncid, dim_name, ncdim_type, ncdim_id )
5633	CALL netcdf_handle_error( 'netcdf_create_dim', error_no )
5634
5635	#endif
5636
5637	END SUBROUTINE netcdf_create_dim
5638
5639
5640	!------------------------------------------------------------------------------!
5641	! Description:
5642	! ------------
5643	!> Create a one dimensional variable in specific units in NetCDF file
5644	!------------------------------------------------------------------------------!
5645
5646	SUBROUTINE netcdf_create_var( ncid, dim_id, var_name, var_type, var_id, &
5647	unit_name, long_name, error_no1, error_no2, &
5648	error_no3, fill )
5649
5650	#if defined( __netcdf )
5651	IMPLICIT NONE
5652
5653	CHARACTER(LEN=*), INTENT(IN) :: long_name
5654	CHARACTER(LEN=*), INTENT(IN) :: unit_name
5655	CHARACTER(LEN=*), INTENT(IN) :: var_name
5656
5657	LOGICAL, OPTIONAL :: fill !< indicates setting of _FillValue attribute
5658
5659	INTEGER, INTENT(IN) :: error_no1
5660	INTEGER, INTENT(IN) :: error_no2
5661	INTEGER, INTENT(IN) :: error_no3
5662	INTEGER, INTENT(IN) :: ncid
5663	INTEGER, INTENT(OUT) :: var_id
5664	INTEGER, INTENT(IN) :: var_type
5665
5666	INTEGER, DIMENSION(:), INTENT(IN) :: dim_id
5667
5668	!
5669	!-- Define variable
5670	nc_stat = NF90_DEF_VAR( ncid, var_name, var_type, dim_id, var_id )
5671	CALL netcdf_handle_error( 'netcdf_create_var', error_no1 )
5672
5673	#if defined( __netcdf4 )
5674	!
5675	!-- Check if variable should be deflate (including shuffling)
5676	!-- and if it is possible (only NetCDF4 with HDF5 supports compression)
5677	IF ( netcdf_data_format > 2 .AND. netcdf_deflate > 0 ) THEN
5678	nc_stat = NF90_DEF_VAR_DEFLATE( ncid, var_id, 1, 1, netcdf_deflate )
5679	CALL netcdf_handle_error( 'netcdf_create_var_deflate', error_no1 )
5680	ENDIF
5681	#endif
5682	!
5683	!-- Set unit name if set
5684	IF ( unit_name /= '' ) THEN
5685	nc_stat = NF90_PUT_ATT( ncid, var_id, 'units', unit_name )
5686	CALL netcdf_handle_error( 'netcdf_create_var', error_no2 )
5687	ENDIF
5688
5689	!
5690	!-- Set long name if set
5691	IF ( long_name /= '' ) THEN
5692	nc_stat = NF90_PUT_ATT( ncid, var_id, 'long_name', long_name )
5693	CALL netcdf_handle_error( 'netcdf_create_var', error_no3 )
5694	ENDIF
5695
5696	!
5697	!-- Set _FillValue for all variables, except for dimension variables.
5698	!-- Set the fill values accordingly to the corresponding output precision.
5699	IF ( PRESENT( fill ) ) THEN
5700	IF ( var_type == NF90_REAL4 ) THEN
5701	nc_stat = NF90_PUT_ATT( ncid, var_id, '_FillValue', &
5702	REAL( fill_value, KIND = 4 ) )
5703	CALL netcdf_handle_error( 'netcdf_create_var', 0 )
5704	ELSE
5705	nc_stat = NF90_PUT_ATT( ncid, var_id, '_FillValue', &
5706	REAL( fill_value, KIND = 8 ) )
5707	CALL netcdf_handle_error( 'netcdf_create_var', 0 )
5708	ENDIF
5709	ENDIF
5710
5711	#endif
5712	END SUBROUTINE netcdf_create_var
5713
5714	END MODULE netcdf_interface

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |