Home

Context Navigation

netcdf_interface_mod.f90 @ 1976

Last change on this file since 1976 was 1976, checked in by maronga, 7 years ago
further modularization of land surface model (2D/3D output and restart data). Bugfix for restart runs without land surface model
Property svn:keywords set to `Id`
File size: 230.8 KB

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