Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |