Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |