Home

Context Navigation

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

Last change on this file since 220 was 216, checked in by raasch, 15 years ago
reading mechanism for restart files completely revised
Property svn:keywords set to `Id`
File size: 143.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |