Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |