Home

Context Navigation

netcdf_interface.f90 @ 1783

Last change on this file since 1783 was 1783, checked in by raasch, 7 years ago
NetCDF routines modularized; new parameter netcdf_deflate; further changes in the pmc
Property svn:keywords set to `Id`
File size: 218.7 KB

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