Home

Context Navigation

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

Last change on this file since 205 was 198, checked in by raasch, 16 years ago
file headers updated for the next release 3.5
Property svn:keywords set to `Id`
File size: 143.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |