Home

Context Navigation

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

Last change on this file since 389 was 359, checked in by letzel, 15 years ago
Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if updated attributes are larger than their original size, NF90_PUT_ATT is called in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible performance loss; an alternative strategy would be to ensure equal attribute size in a job chain. (netcdf) Changed: For extended NetCDF files, the updated title attribute includes an update of time_average_text where appropriate. (netcdf) Bugfix: correction of initial volume flow for non-flat topography (init_3d_model) Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model) New: MPI debug option (mrun) Documentation updates: topography file format (chapter 3.5.4) NC identifier range expanded to 444 (appendix_a)
Property svn:keywords set to `Id`
File size: 149.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |