Home

Context Navigation

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

Last change on this file since 138 was 98, checked in by raasch, 18 years ago
updating comments and rc-file
Property svn:keywords set to `Id`
File size: 142.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |