Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |