Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |