Home

Context Navigation

netcdf.f90 @ 1763

Last change on this file since 1763 was 1746, checked in by gronemeier, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 219.2 KB

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

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

Context Navigation

source: palm/trunk/SOURCE/netcdf.f90 @ 1763

Download in other formats: