Home

Context Navigation

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

Last change on this file since 1025 was 993, checked in by hoffmann, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 184.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |