Home

Context Navigation

netcdf.f90 @ 145

Last change on this file since 145 was 144, checked in by letzel, 17 years ago

User-defined spectra.

Bugfix: extra '*' removed in user_statistics sample code (user_interface).

Property svn:keywords set to Id

File size: 143.5 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	! user-defined spectra
10	!
11	! Former revisions:
12	! -----------------
13	! $Id: netcdf.f90 144 2008-01-04 04:29:45Z raasch $
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	!-- First check for the vertical grid
2871	found = .TRUE.
2872	SELECT CASE ( data_output_sp(i) )
2873	!
2874	!-- Most variables are defined on the zu levels
2875	CASE ( 'e', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c', 'ql_v', &
2876	'ql_vp', 'qv', 'rho', 's', 'sa', 'u', 'v', 'vpt' )
2877
2878	grid_z = 'zu'
2879	!
2880	!-- zw levels
2881	CASE ( 'w' )
2882
2883	grid_z = 'zw'
2884
2885	CASE DEFAULT
2886	!
2887	!-- Check for user-defined quantities (found, grid_x and grid_y
2888	!-- are dummies)
2889	CALL user_define_netcdf_grid( data_output_sp(i), found, &
2890	grid_x, grid_y, grid_z )
2891
2892	END SELECT
2893
2894	IF ( INDEX( spectra_direction(i), 'x' ) /= 0 ) THEN
2895
2896	!
2897	!-- Define the variable
2898	netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
2899	CALL clean_netcdf_varname( netcdf_var_name )
2900	IF ( TRIM( grid_z ) == 'zw' ) THEN
2901	nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name, &
2902	nc_precision(7), (/ id_dim_x_sp, &
2903	id_dim_zw_sp, id_dim_time_sp /), &
2904	id_var_dospx(i) )
2905	ELSE
2906	nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name, &
2907	nc_precision(7), (/ id_dim_x_sp, &
2908	id_dim_zu_sp, id_dim_time_sp /), &
2909	id_var_dospx(i) )
2910	ENDIF
2911	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 285 )
2912
2913	nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospx(i), &
2914	'long_name', netcdf_var_name )
2915	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 286 )
2916
2917	nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospx(i), &
2918	'units', 'unknown' )
2919	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 287 )
2920
2921	var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
2922
2923	ENDIF
2924
2925	IF ( INDEX( spectra_direction(i), 'y' ) /= 0 ) THEN
2926
2927	!
2928	!-- Define the variable
2929	netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
2930	CALL clean_netcdf_varname( netcdf_var_name )
2931	IF ( TRIM( grid_z ) == 'zw' ) THEN
2932	nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name, &
2933	nc_precision(7), (/ id_dim_y_sp, &
2934	id_dim_zw_sp, id_dim_time_sp /), &
2935	id_var_dospy(i) )
2936	ELSE
2937	nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name, &
2938	nc_precision(7), (/ id_dim_y_sp, &
2939	id_dim_zu_sp, id_dim_time_sp /), &
2940	id_var_dospy(i) )
2941	ENDIF
2942	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 288 )
2943
2944	nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospy(i), &
2945	'long_name', netcdf_var_name )
2946	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 289 )
2947
2948	nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospy(i), &
2949	'units', 'unknown' )
2950	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 290 )
2951
2952	var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
2953
2954	ENDIF
2955
2956	i = i + 1
2957
2958	ENDDO
2959
2960	!
2961	!-- Write the list of variables as global attribute (this is used by
2962	!-- restart runs)
2963	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', var_list )
2964	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 291 )
2965
2966	!
2967	!-- Leave NetCDF define mode
2968	nc_stat = NF90_ENDDEF( id_set_sp )
2969	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 292 )
2970
2971	!
2972	!-- Write axis data: zu_sp, zw_sp, k_x, k_y
2973	ALLOCATE( netcdf_data(1:ns) )
2974
2975	!
2976	!-- Write zu data
2977	netcdf_data(1:ns) = zu( comp_spectra_level(1:ns) )
2978	nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zu_sp, netcdf_data, &
2979	start = (/ 1 /), count = (/ ns /) )
2980	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 293 )
2981
2982	!
2983	!-- Write zw data
2984	netcdf_data(1:ns) = zw( comp_spectra_level(1:ns) )
2985	nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zw_sp, netcdf_data, &
2986	start = (/ 1 /), count = (/ ns /) )
2987	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 294 )
2988
2989	DEALLOCATE( netcdf_data )
2990
2991	!
2992	!-- Write data for x and y axis (wavenumbers)
2993	ALLOCATE( netcdf_data(nx/2) )
2994	DO i = 1, nx/2
2995	netcdf_data(i) = 2.0 * pi * i / ( dx * ( nx + 1 ) )
2996	ENDDO
2997
2998	nc_stat = NF90_PUT_VAR( id_set_sp, id_var_x_sp, netcdf_data, &
2999	start = (/ 1 /), count = (/ nx/2 /) )
3000	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 295 )
3001
3002	DEALLOCATE( netcdf_data )
3003
3004	ALLOCATE( netcdf_data(ny/2) )
3005	DO i = 1, ny/2
3006	netcdf_data(i) = 2.0 * pi * i / ( dy * ( ny + 1 ) )
3007	ENDDO
3008
3009	nc_stat = NF90_PUT_VAR( id_set_sp, id_var_y_sp, netcdf_data, &
3010	start = (/ 1 /), count = (/ ny/2 /) )
3011	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 296 )
3012
3013	DEALLOCATE( netcdf_data )
3014
3015
3016	CASE ( 'sp_ext' )
3017
3018	!
3019	!-- Get the list of variables and compare with the actual run.
3020	!-- First var_list_old has to be reset, since GET_ATT does not assign
3021	!-- trailing blanks.
3022	var_list_old = ' '
3023	nc_stat = NF90_GET_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', &
3024	var_list_old )
3025	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 297 )
3026
3027	var_list = ';'
3028	i = 1
3029	DO WHILE ( data_output_sp(i) /= ' ' .AND. i <= 10 )
3030
3031	IF ( INDEX( spectra_direction(i), 'x' ) /= 0 ) THEN
3032	netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
3033	var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
3034	ENDIF
3035
3036	IF ( INDEX( spectra_direction(i), 'y' ) /= 0 ) THEN
3037	netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
3038	var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
3039	ENDIF
3040
3041	i = i + 1
3042
3043	ENDDO
3044
3045	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
3046	PRINT*, '+++ WARNING: NetCDF file for spectra from previous ' //&
3047	'run found,'
3048	PRINT*, ' but this file cannot be extended due to' // &
3049	' variable mismatch.'
3050	PRINT*, ' New file is created instead.'
3051	extend = .FALSE.
3052	RETURN
3053	ENDIF
3054
3055	!
3056	!-- Determine the number of current vertical levels for which spectra
3057	!-- shall be output
3058	ns = 1
3059	DO WHILE ( comp_spectra_level(ns) /= 999999 .AND. ns <= 100 )
3060	ns = ns + 1
3061	ENDDO
3062	ns = ns - 1
3063
3064	!
3065	!-- Get and compare the number of vertical levels
3066	nc_stat = NF90_INQ_VARID( id_set_sp, 'zu_sp', id_var_zu_sp )
3067	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 298 )
3068
3069	nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_zu_sp, &
3070	dimids = id_dim_zu_sp_old )
3071	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 299 )
3072	id_dim_zu_sp = id_dim_zu_sp_old(1)
3073
3074	nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_zu_sp, &
3075	len = ns_old )
3076	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 300 )
3077
3078	IF ( ns /= ns_old ) THEN
3079	PRINT*, '+++ WARNING: NetCDF file for spectra from previous ' //&
3080	'run found,'
3081	PRINT*, ' but this file cannot be extended due to' // &
3082	' mismatch in number of'
3083	PRINT*, ' vertical levels.'
3084	PRINT*, ' New file is created instead.'
3085	extend = .FALSE.
3086	RETURN
3087	ENDIF
3088
3089	!
3090	!-- Get and compare the heights of the cross sections
3091	ALLOCATE( netcdf_data(1:ns_old) )
3092
3093	nc_stat = NF90_GET_VAR( id_set_sp, id_var_zu_sp, netcdf_data )
3094	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 301 )
3095
3096	DO i = 1, ns
3097	IF ( zu(comp_spectra_level(i)) /= netcdf_data(i) ) THEN
3098	PRINT*, '+++ WARNING: NetCDF file for spectra from previou' // &
3099	's run found,'
3100	PRINT*, ' but this file cannot be extended due' // &
3101	' to mismatch in heights'
3102	PRINT*, ' of vertical levels.'
3103	PRINT*, ' New file is created instead.'
3104	extend = .FALSE.
3105	RETURN
3106	ENDIF
3107	ENDDO
3108
3109	DEALLOCATE( netcdf_data )
3110
3111	!
3112	!-- Get the id of the time coordinate (unlimited coordinate) and its
3113	!-- last index on the file. The next time level is plsp..count+1.
3114	!-- The current time must be larger than the last output time
3115	!-- on the file.
3116	nc_stat = NF90_INQ_VARID( id_set_sp, 'time', id_var_time_sp )
3117	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 302 )
3118
3119	nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_time_sp, &
3120	dimids = id_dim_time_old )
3121	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 303 )
3122	id_dim_time_sp = id_dim_time_old(1)
3123
3124	nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_time_sp, &
3125	len = dosp_time_count )
3126	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 304 )
3127
3128	nc_stat = NF90_GET_VAR( id_set_sp, id_var_time_sp, &
3129	last_time_coordinate, &
3130	start = (/ dosp_time_count /), &
3131	count = (/ 1 /) )
3132	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 305 )
3133
3134	IF ( last_time_coordinate(1) >= simulated_time ) THEN
3135	PRINT*, '+++ WARNING: NetCDF file for spectra from previous ' // &
3136	'run found,'
3137	PRINT*, ' but this file cannot be extended becaus' // &
3138	'e the current output time'
3139	PRINT*, ' is less or equal than the last output t' // &
3140	'ime on this file.'
3141	PRINT*, ' New file is created instead.'
3142	dosp_time_count = 0
3143	extend = .FALSE.
3144	RETURN
3145	ENDIF
3146
3147	!
3148	!-- Dataset seems to be extendable.
3149	!-- Now get the variable ids.
3150	i = 1
3151	DO WHILE ( data_output_sp(i) /= ' ' .AND. i <= 10 )
3152
3153	IF ( INDEX( spectra_direction(i), 'x' ) /= 0 ) THEN
3154	netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
3155	nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
3156	id_var_dospx(i) )
3157	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 306 )
3158	ENDIF
3159
3160	IF ( INDEX( spectra_direction(i), 'y' ) /= 0 ) THEN
3161	netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
3162	nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
3163	id_var_dospy(i) )
3164	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 307 )
3165	ENDIF
3166
3167	i = i + 1
3168
3169	ENDDO
3170
3171	!
3172	!-- Change the titel attribute on file
3173	IF ( averaging_interval_sp /= 0.0 ) THEN
3174	WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
3175	averaging_interval_sp
3176	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
3177	TRIM( run_description_header ) // &
3178	TRIM( time_average_text ) )
3179	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 308 )
3180
3181	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval_sp
3182	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'time_avg', &
3183	TRIM( time_average_text ) )
3184	ELSE
3185	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
3186	TRIM( run_description_header ) )
3187	ENDIF
3188	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 309 )
3189
3190	PRINT, '** NetCDF file for spectra from previous run found.'
3191	PRINT*, ' This file will be extended.'
3192
3193
3194	CASE ( 'pt_new' )
3195
3196	!
3197	!-- Define some global attributes of the dataset
3198	nc_stat = NF90_PUT_ATT( id_set_prt, NF90_GLOBAL, 'title', &
3199	TRIM( run_description_header ) )
3200	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 310 )
3201
3202	!
3203	!-- Define time coordinate for particles (unlimited dimension)
3204	nc_stat = NF90_DEF_DIM( id_set_prt, 'time', NF90_UNLIMITED, &
3205	id_dim_time_prt )
3206	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 311 )
3207
3208	nc_stat = NF90_DEF_VAR( id_set_prt, 'time', NF90_DOUBLE, &
3209	id_dim_time_prt, id_var_time_prt )
3210	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 312 )
3211
3212	nc_stat = NF90_PUT_ATT( id_set_prt, id_var_time_prt, 'units', &
3213	'seconds' )
3214	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 313 )
3215
3216	!
3217	!-- Define particle coordinate (maximum particle number)
3218	nc_stat = NF90_DEF_DIM( id_set_prt, 'particle_number', &
3219	maximum_number_of_particles, id_dim_prtnum )
3220	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 314 )
3221
3222	nc_stat = NF90_DEF_VAR( id_set_prt, 'particle_number', NF90_DOUBLE, &
3223	id_dim_prtnum, id_var_prtnum )
3224	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 315 )
3225
3226	nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prtnum, 'units', &
3227	'particle number' )
3228	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 316 )
3229
3230	!
3231	!-- Define variable which contains the real number of particles in use
3232	nc_stat = NF90_DEF_VAR( id_set_prt, 'real_num_of_prt', NF90_INT, &
3233	id_dim_time_prt, id_var_rnop_prt )
3234	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 317 )
3235
3236	nc_stat = NF90_PUT_ATT( id_set_prt, id_var_rnop_prt, 'units', &
3237	'particle number' )
3238	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 318 )
3239
3240	!
3241	!-- Define the variables
3242	DO i = 1, 17
3243
3244	nc_stat = NF90_DEF_VAR( id_set_prt, prt_var_names(i), &
3245	nc_precision(8), &
3246	(/ id_dim_prtnum, id_dim_time_prt /), &
3247	id_var_prt(i) )
3248	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 319 )
3249
3250	nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prt(i), &
3251	'long_name', TRIM( prt_var_names(i) ) )
3252	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 320 )
3253
3254	nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prt(i), &
3255	'units', TRIM( prt_var_units(i) ) )
3256	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 321 )
3257
3258	ENDDO
3259
3260	!
3261	!-- Leave NetCDF define mode
3262	nc_stat = NF90_ENDDEF( id_set_prt )
3263	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 322 )
3264
3265
3266	CASE ( 'pt_ext' )
3267
3268	!
3269	!-- Get the id of the time coordinate (unlimited coordinate) and its
3270	!-- last index on the file. The next time level is prt..count+1.
3271	!-- The current time must be larger than the last output time
3272	!-- on the file.
3273	nc_stat = NF90_INQ_VARID( id_set_prt, 'time', id_var_time_prt )
3274	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 323 )
3275
3276	nc_stat = NF90_INQUIRE_VARIABLE( id_set_prt, id_var_time_prt, &
3277	dimids = id_dim_time_old )
3278	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 324 )
3279	id_dim_time_prt = id_dim_time_old(1)
3280
3281	nc_stat = NF90_INQUIRE_DIMENSION( id_set_prt, id_dim_time_prt, &
3282	len = prt_time_count )
3283	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 325 )
3284
3285	nc_stat = NF90_GET_VAR( id_set_prt, id_var_time_prt, &
3286	last_time_coordinate, &
3287	start = (/ prt_time_count /), &
3288	count = (/ 1 /) )
3289	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 326 )
3290
3291	IF ( last_time_coordinate(1) >= simulated_time ) THEN
3292	PRINT*, '+++ WARNING: NetCDF file for particles from previous ' //&
3293	'run found,'
3294	PRINT*, ' but this file cannot be extended becaus' // &
3295	'e the current output time'
3296	PRINT*, ' is less or equal than the last output t' // &
3297	'ime on this file.'
3298	PRINT*, ' New file is created instead.'
3299	prt_time_count = 0
3300	extend = .FALSE.
3301	RETURN
3302	ENDIF
3303
3304	!
3305	!-- Dataset seems to be extendable.
3306	!-- Now get the variable ids.
3307	nc_stat = NF90_INQ_VARID( id_set_prt, 'real_num_of_prt', &
3308	id_var_rnop_prt )
3309	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 327 )
3310
3311	DO i = 1, 17
3312
3313	nc_stat = NF90_INQ_VARID( id_set_prt, prt_var_names(i), &
3314	id_var_prt(i) )
3315	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 328 )
3316
3317	ENDDO
3318
3319	IF ( myid == 0 ) THEN
3320	PRINT, '** NetCDF file for particles from previous run found.'
3321	PRINT*, ' This file will be extended.'
3322	ENDIF
3323
3324
3325	CASE ( 'ps_new' )
3326
3327	!
3328	!-- Define some global attributes of the dataset
3329	nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
3330	TRIM( run_description_header ) )
3331	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 396 )
3332
3333	!
3334	!-- Define time coordinate for particle time series (unlimited dimension)
3335	nc_stat = NF90_DEF_DIM( id_set_pts, 'time', NF90_UNLIMITED, &
3336	id_dim_time_pts )
3337	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 397 )
3338
3339	nc_stat = NF90_DEF_VAR( id_set_pts, 'time', NF90_DOUBLE, &
3340	id_dim_time_pts, id_var_time_pts )
3341	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 398 )
3342
3343	nc_stat = NF90_PUT_ATT( id_set_pts, id_var_time_pts, 'units', &
3344	'seconds')
3345	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 399 )
3346
3347	!
3348	!-- Define the variables. If more than one particle group is defined,
3349	!-- define seperate variables for each group
3350	var_list = ';'
3351	DO i = 1, dopts_num
3352
3353	!
3354	!-- First, remove those characters not allowed by NetCDF
3355	netcdf_var_name = dopts_label(i)
3356	CALL clean_netcdf_varname( netcdf_var_name )
3357
3358	DO j = 0, number_of_particle_groups
3359
3360	IF ( j == 0 ) THEN
3361	suffix1 = ''
3362	ELSE
3363	WRITE ( suffix1, '(''_'',I2.2)' ) j
3364	ENDIF
3365
3366	nc_stat = NF90_DEF_VAR( id_set_pts, &
3367	TRIM( netcdf_var_name ) // suffix1, &
3368	nc_precision(6), &
3369	(/ id_dim_time_pts /), &
3370	id_var_dopts(i,j) )
3371	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 400 )
3372
3373	IF ( j == 0 ) THEN
3374	nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
3375	'long_name', &
3376	TRIM( dopts_label(i) ) )
3377	ELSE
3378	nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
3379	'long_name', &
3380	TRIM( dopts_label(i) ) // ' PG ' // &
3381	suffix1(2:3) )
3382	ENDIF
3383	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 401 )
3384
3385	nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
3386	'units', TRIM( dopts_unit(i) ) )
3387	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 402 )
3388
3389	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
3390	suffix1 // ';'
3391
3392	IF ( number_of_particle_groups == 1 ) EXIT
3393
3394	ENDDO
3395
3396	ENDDO
3397
3398	!
3399	!-- Write the list of variables as global attribute (this is used by
3400	!-- restart runs)
3401	nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
3402	var_list )
3403	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 403 )
3404
3405	!
3406	!-- Leave NetCDF define mode
3407	nc_stat = NF90_ENDDEF( id_set_pts )
3408	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 404 )
3409
3410
3411	CASE ( 'ps_ext' )
3412
3413	!
3414	!-- Get the list of variables and compare with the actual run.
3415	!-- First var_list_old has to be reset, since GET_ATT does not assign
3416	!-- trailing blanks.
3417	var_list_old = ' '
3418	nc_stat = NF90_GET_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
3419	var_list_old )
3420	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 405 )
3421
3422	var_list = ';'
3423	i = 1
3424	DO i = 1, dopts_num
3425
3426	netcdf_var_name = dopts_label(i)
3427	CALL clean_netcdf_varname( netcdf_var_name )
3428
3429	DO j = 0, number_of_particle_groups
3430
3431	IF ( j == 0 ) THEN
3432	suffix1 = ''
3433	ELSE
3434	WRITE ( suffix1, '(''_'',I2.2)' ) j
3435	ENDIF
3436
3437	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
3438	suffix1 // ';'
3439
3440	IF ( number_of_particle_groups == 1 ) EXIT
3441
3442	ENDDO
3443
3444	ENDDO
3445
3446	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
3447	PRINT*, '+++ WARNING: NetCDF file for particle time series ' //&
3448	'from previuos run found,'
3449	PRINT*, ' but this file cannot be extended due to' // &
3450	' variable mismatch.'
3451	PRINT*, ' New file is created instead.'
3452	extend = .FALSE.
3453	RETURN
3454	ENDIF
3455
3456	!
3457	!-- Get the id of the time coordinate (unlimited coordinate) and its
3458	!-- last index on the file. The next time level is dots..count+1.
3459	!-- The current time must be larger than the last output time
3460	!-- on the file.
3461	nc_stat = NF90_INQ_VARID( id_set_pts, 'time', id_var_time_pts )
3462	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 406 )
3463
3464	nc_stat = NF90_INQUIRE_VARIABLE( id_set_pts, id_var_time_pts, &
3465	dimids = id_dim_time_old )
3466	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 407 )
3467	id_dim_time_pts = id_dim_time_old(1)
3468
3469	nc_stat = NF90_INQUIRE_DIMENSION( id_set_pts, id_dim_time_pts, &
3470	len = dopts_time_count )
3471	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 408 )
3472
3473	nc_stat = NF90_GET_VAR( id_set_pts, id_var_time_pts, &
3474	last_time_coordinate, &
3475	start = (/ dopts_time_count /), &
3476	count = (/ 1 /) )
3477	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 409 )
3478
3479	IF ( last_time_coordinate(1) >= simulated_time ) THEN
3480	PRINT*, '+++ WARNING: NetCDF file for time series from' // &
3481	' previuos run found,'
3482	PRINT*, ' but this file cannot be extended becaus' // &
3483	'e the current output time'
3484	PRINT*, ' is less or equal than the last output t' // &
3485	'ime on this file.'
3486	PRINT*, ' New file is created instead.'
3487	dopts_time_count = 0
3488	extend = .FALSE.
3489	RETURN
3490	ENDIF
3491
3492	!
3493	!-- Dataset seems to be extendable.
3494	!-- Now get the variable ids
3495	i = 1
3496	DO i = 1, dopts_num
3497
3498	netcdf_var_name_base = dopts_label(i)
3499	CALL clean_netcdf_varname( netcdf_var_name_base )
3500
3501	DO j = 0, number_of_particle_groups
3502
3503	IF ( j == 0 ) THEN
3504	suffix1 = ''
3505	ELSE
3506	WRITE ( suffix1, '(''_'',I2.2)' ) j
3507	ENDIF
3508
3509	netcdf_var_name = TRIM( netcdf_var_name_base ) // suffix1
3510
3511	nc_stat = NF90_INQ_VARID( id_set_pts, netcdf_var_name, &
3512	id_var_dopts(i,j) )
3513	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 410 )
3514
3515	IF ( number_of_particle_groups == 1 ) EXIT
3516
3517	ENDDO
3518
3519	ENDDO
3520
3521	!
3522	!-- Change the title attribute on file
3523	nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
3524	TRIM( run_description_header ) )
3525	IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 411 )
3526	PRINT, '** NetCDF file for particle time series from previous ', &
3527	'run found.'
3528	PRINT*, ' This file will be extended.'
3529
3530
3531	CASE DEFAULT
3532
3533	PRINT*, '+++ define_netcdf_header: mode "', mode, '" not supported'
3534
3535	END SELECT
3536
3537	#endif
3538	END SUBROUTINE define_netcdf_header
3539
3540
3541
3542	SUBROUTINE handle_netcdf_error( errno )
3543	#if defined( __netcdf )
3544
3545	!------------------------------------------------------------------------------!
3546	!
3547	! Description:
3548	! ------------
3549	! Prints out a text message corresponding to the current status.
3550	!------------------------------------------------------------------------------!
3551
3552	USE netcdf
3553	USE netcdf_control
3554	USE pegrid
3555
3556	IMPLICIT NONE
3557
3558	INTEGER :: errno
3559
3560	IF ( nc_stat /= NF90_NOERR ) THEN
3561	PRINT*, '+++ netcdf error ', errno,': ', TRIM( NF90_STRERROR( nc_stat ) )
3562	#if defined( __parallel )
3563	CALL MPI_ABORT( comm2d, 9999, ierr )
3564	#else
3565	CALL local_stop
3566	#endif
3567	ENDIF
3568
3569	#endif
3570	END SUBROUTINE handle_netcdf_error
3571
3572
3573
3574	SUBROUTINE clean_netcdf_varname( string )
3575	#if defined( __netcdf )
3576
3577	!------------------------------------------------------------------------------!
3578	!
3579	! Description:
3580	! ------------
3581	! Replace those characters in string which are not allowed by NetCDF.
3582	!------------------------------------------------------------------------------!
3583
3584	USE netcdf_control
3585
3586	IMPLICIT NONE
3587
3588	CHARACTER (LEN=10), INTENT(INOUT) :: string
3589
3590	INTEGER :: i, ic
3591
3592	DO i = 1, replace_num
3593	DO
3594	ic = INDEX( string, replace_char(i) )
3595	IF ( ic /= 0 ) THEN
3596	string(ic:ic) = replace_by(i)
3597	ELSE
3598	EXIT
3599	ENDIF
3600	ENDDO
3601	ENDDO
3602
3603	#endif
3604	END SUBROUTINE clean_netcdf_varname

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

Context Navigation

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

Download in other formats: