Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |