Home

Context Navigation

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

Last change on this file since 510 was 494, checked in by raasch, 14 years ago
last commit documented; configuration example file for netcdf4 added
Property svn:keywords set to `Id`
File size: 182.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 494 2010-03-01 12:23:32Z raasch $
14	!
15	! 493 2010-03-01 08:30:24Z raasch
16	! Extensions for NetCDF4 output
17	!
18	! 410 2009-12-04 17:05:40Z letzel
19	! masked data output
20	!
21	! 359 2009-08-19 16:56:44Z letzel
22	! For extended NetCDF files, the updated title attribute includes an update of
23	! time_average_text where appropriate.
24	! Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
25	! updated attributes are larger than their original size, NF90_PUT_ATT is called
26	! in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
27	! possible performance loss; an alternative strategy would be to ensure equal
28	! attribute size in a job chain.
29	! NetCDF unit attribute in timeseries output in case of statistic
30	! regions added.
31	! Output of NetCDF messages with aid of message handling routine.
32	! Output of messages replaced by message handling routine.
33	! Typographical errors fixed.
34	!
35	! 216 2008-11-25 07:12:43Z raasch
36	! Origin of the xy-coordinate system shifted from the center of the first
37	! grid cell (indices i=0, j=0) to the south-left corner of this cell.
38	!
39	! 189 2008-08-13 17:09:26Z letzel
40	! consistently allow 100 spectra levels instead of 10
41	! bug fix in the determination of the number of output heights for spectra,
42	! +user-defined spectra
43	!
44	! 97 2007-06-21 08:23:15Z raasch
45	! Grids defined for rho and sa
46	!
47	! 48 2007-03-06 12:28:36Z raasch
48	! Output topography height information (zu_s_inner, zw_s_inner) to 2d-xy and 3d
49	! datasets
50	!
51	! RCS Log replace by Id keyword, revision history cleaned up
52	!
53	! Revision 1.12 2006/09/26 19:35:16 raasch
54	! Bugfix yv coordinates for yz cross sections
55	!
56	! Revision 1.1 2005/05/18 15:37:16 raasch
57	! Initial revision
58	!
59	!
60	! Description:
61	! ------------
62	! In case of extend = .FALSE.:
63	! Define all necessary dimensions, axes and variables for the different
64	! NetCDF datasets. This subroutine is called from check_open after a new
65	! dataset is created. It leaves the open NetCDF files ready to write.
66	!
67	! In case of extend = .TRUE.:
68	! Find out if dimensions and variables of an existing file match the values
69	! of the actual run. If so, get all necessary informations (ids, etc.) from
70	! this file.
71	!
72	! Parameter av can assume values 0 (non-averaged data) and 1 (time averaged
73	! data)
74	!------------------------------------------------------------------------------!
75	#if defined( __netcdf )
76
77	USE arrays_3d
78	USE constants
79	USE control_parameters
80	USE grid_variables
81	USE indices
82	USE netcdf_control
83	USE pegrid
84	USE particle_attributes
85	USE profil_parameter
86	USE spectrum
87	USE statistics
88
89
90	IMPLICIT NONE
91
92	CHARACTER (LEN=2) :: suffix
93	CHARACTER (LEN=2), INTENT (IN) :: callmode
94	CHARACTER (LEN=3) :: suffix1
95	CHARACTER (LEN=4) :: grid_x, grid_y, grid_z
96	CHARACTER (LEN=6) :: mode
97	CHARACTER (LEN=10) :: netcdf_var_name, netcdf_var_name_base, &
98	precision, var
99	CHARACTER (LEN=80) :: time_average_text
100	CHARACTER (LEN=2000) :: var_list, var_list_old
101
102	INTEGER :: av, file_id, i, id_x, id_y, id_z, j, ns, ns_old, nz_old
103
104	INTEGER, DIMENSION(1) :: id_dim_time_old, id_dim_x_yz_old, &
105	id_dim_y_xz_old, id_dim_zu_sp_old, &
106	id_dim_zu_xy_old, id_dim_zu_3d_old, &
107	id_dim_zu_mask_old
108
109	LOGICAL :: found
110
111	LOGICAL, INTENT (INOUT) :: extend
112
113	LOGICAL, SAVE :: init_netcdf = .FALSE.
114
115	REAL, DIMENSION(1) :: last_time_coordinate
116
117	REAL, DIMENSION(:), ALLOCATABLE :: netcdf_data
118	REAL, DIMENSION(:,:), ALLOCATABLE :: netcdf_data_2d
119
120
121
122	!
123	!-- Initializing actions (return to calling routine check_parameters afterwards)
124	IF ( .NOT. init_netcdf ) THEN
125	!
126	!-- Check and set accuracy for NetCDF output. First set default value
127	nc_precision = NF90_REAL4
128
129	i = 1
130	DO WHILE ( netcdf_precision(i) /= ' ' )
131	j = INDEX( netcdf_precision(i), '_' )
132	IF ( j == 0 ) THEN
133	WRITE ( message_string, * ) 'netcdf_precision must contain a ', &
134	'"_"netcdf_precision(', i, ')="', &
135	TRIM( netcdf_precision(i) ),'"'
136	CALL message( 'define_netcdf_header', 'PA0241', 1, 2, 0, 6, 0 )
137	ENDIF
138
139	var = netcdf_precision(i)(1:j-1)
140	precision = netcdf_precision(i)(j+1:)
141
142	IF ( precision == 'NF90_REAL4' ) THEN
143	j = NF90_REAL4
144	ELSEIF ( precision == 'NF90_REAL8' ) THEN
145	j = NF90_REAL8
146	ELSE
147	WRITE ( message_string, * ) 'illegal netcdf precision: ', &
148	'netcdf_precision(', i, ')="', &
149	TRIM( netcdf_precision(i) ),'"'
150	CALL message( 'define_netcdf_header', 'PA0242', 1, 2, 0, 6, 0 )
151	ENDIF
152
153	SELECT CASE ( var )
154	CASE ( 'xy' )
155	nc_precision(1) = j
156	CASE ( 'xz' )
157	nc_precision(2) = j
158	CASE ( 'yz' )
159	nc_precision(3) = j
160	CASE ( '2d' )
161	nc_precision(1:3) = j
162	CASE ( '3d' )
163	nc_precision(4) = j
164	CASE ( 'pr' )
165	nc_precision(5) = j
166	CASE ( 'ts' )
167	nc_precision(6) = j
168	CASE ( 'sp' )
169	nc_precision(7) = j
170	CASE ( 'prt' )
171	nc_precision(8) = j
172	CASE ( 'masks' )
173	nc_precision(11:50) = j
174	CASE ( 'mask01' )
175	nc_precision(11) = j
176	CASE ( 'mask02' )
177	nc_precision(12) = j
178	CASE ( 'mask03' )
179	nc_precision(13) = j
180	CASE ( 'mask04' )
181	nc_precision(14) = j
182	CASE ( 'mask05' )
183	nc_precision(15) = j
184	CASE ( 'mask06' )
185	nc_precision(16) = j
186	CASE ( 'mask07' )
187	nc_precision(17) = j
188	CASE ( 'mask08' )
189	nc_precision(18) = j
190	CASE ( 'mask09' )
191	nc_precision(19) = j
192	CASE ( 'mask10' )
193	nc_precision(20) = j
194	CASE ( 'mask11' )
195	nc_precision(21) = j
196	CASE ( 'mask12' )
197	nc_precision(22) = j
198	CASE ( 'mask13' )
199	nc_precision(23) = j
200	CASE ( 'mask14' )
201	nc_precision(24) = j
202	CASE ( 'mask15' )
203	nc_precision(25) = j
204	CASE ( 'mask16' )
205	nc_precision(26) = j
206	CASE ( 'mask17' )
207	nc_precision(27) = j
208	CASE ( 'mask18' )
209	nc_precision(28) = j
210	CASE ( 'mask19' )
211	nc_precision(29) = j
212	CASE ( 'mask20' )
213	nc_precision(30) = j
214	CASE ( 'maskav01' )
215	nc_precision(31) = j
216	CASE ( 'maskav02' )
217	nc_precision(32) = j
218	CASE ( 'maskav03' )
219	nc_precision(33) = j
220	CASE ( 'maskav04' )
221	nc_precision(34) = j
222	CASE ( 'maskav05' )
223	nc_precision(35) = j
224	CASE ( 'maskav06' )
225	nc_precision(36) = j
226	CASE ( 'maskav07' )
227	nc_precision(37) = j
228	CASE ( 'maskav08' )
229	nc_precision(38) = j
230	CASE ( 'maskav09' )
231	nc_precision(39) = j
232	CASE ( 'maskav10' )
233	nc_precision(40) = j
234	CASE ( 'maskav11' )
235	nc_precision(41) = j
236	CASE ( 'maskav12' )
237	nc_precision(42) = j
238	CASE ( 'maskav13' )
239	nc_precision(43) = j
240	CASE ( 'maskav14' )
241	nc_precision(44) = j
242	CASE ( 'maskav15' )
243	nc_precision(45) = j
244	CASE ( 'maskav16' )
245	nc_precision(46) = j
246	CASE ( 'maskav17' )
247	nc_precision(47) = j
248	CASE ( 'maskav18' )
249	nc_precision(48) = j
250	CASE ( 'maskav19' )
251	nc_precision(49) = j
252	CASE ( 'maskav20' )
253	nc_precision(50) = j
254	CASE ( 'all' )
255	nc_precision = j
256
257	CASE DEFAULT
258	WRITE ( message_string, * ) 'unknown variable in inipar ', &
259	'assignment: netcdf_precision(', i, ')="', &
260	TRIM( netcdf_precision(i) ),'"'
261	CALL message( 'define_netcdf_header', 'PA0243', 1, 2, 0, 6, 0 )
262
263	END SELECT
264
265	i = i + 1
266	IF ( i > 50 ) EXIT
267	ENDDO
268
269	init_netcdf = .TRUE.
270
271	RETURN
272
273	ENDIF
274
275	!
276	!-- Determine the mode to be processed
277	IF ( extend ) THEN
278	mode = callmode // '_ext'
279	ELSE
280	mode = callmode // '_new'
281	ENDIF
282
283	!
284	!-- Select the mode to be processed. Possibilities are 3d, mask, xy, xz, yz,
285	!-- pr and ts.
286	SELECT CASE ( mode )
287
288	CASE ( 'ma_new' )
289
290	!
291	!-- decompose actual parameter file_id (=formal parameter av) into
292	!-- mid and av
293	file_id = av
294	IF ( file_id <= 140 ) THEN
295	mid = file_id - 120
296	av = 0
297	ELSE
298	mid = file_id - 140
299	av = 1
300	ENDIF
301
302	!
303	!-- Define some global attributes of the dataset
304	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
305	'Conventions', 'COARDS' )
306	CALL handle_netcdf_error( 'netcdf', 9998 )
307
308	IF ( av == 0 ) THEN
309	time_average_text = ' '
310	ELSE
311	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
312	averaging_interval
313	ENDIF
314	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'title', &
315	TRIM( run_description_header ) // &
316	TRIM( time_average_text ) )
317	CALL handle_netcdf_error( 'netcdf', 9998 )
318	IF ( av == 1 ) THEN
319	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval
320	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
321	'time_avg', TRIM( time_average_text ) )
322	CALL handle_netcdf_error( 'netcdf', 9998 )
323	ENDIF
324
325	!
326	!-- Define time coordinate for volume data (unlimited dimension)
327	nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'time', NF90_UNLIMITED, &
328	id_dim_time_mask(mid,av) )
329	CALL handle_netcdf_error( 'netcdf', 9998 )
330
331	nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'time', NF90_DOUBLE, &
332	id_dim_time_mask(mid,av), &
333	id_var_time_mask(mid,av) )
334	CALL handle_netcdf_error( 'netcdf', 9998 )
335
336	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
337	id_var_time_mask(mid,av), 'units', &
338	'seconds')
339	CALL handle_netcdf_error( 'netcdf', 9998 )
340
341	!
342	!-- Define spatial dimensions and coordinates:
343	!-- Define vertical coordinate grid (zu grid)
344	nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'zu_3d', &
345	mask_size(mid,3), id_dim_zu_mask(mid,av) )
346	CALL handle_netcdf_error( 'netcdf', 9998 )
347
348	nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zu_3d', NF90_DOUBLE, &
349	id_dim_zu_mask(mid,av), &
350	id_var_zu_mask(mid,av) )
351	CALL handle_netcdf_error( 'netcdf', 9998 )
352
353	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
354	'units', 'meters' )
355	CALL handle_netcdf_error( 'netcdf', 9998 )
356
357	!
358	!-- Define vertical coordinate grid (zw grid)
359	nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'zw_3d', &
360	mask_size(mid,3), id_dim_zw_mask(mid,av) )
361	CALL handle_netcdf_error( 'netcdf', 9998 )
362
363	nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zw_3d', NF90_DOUBLE, &
364	id_dim_zw_mask(mid,av), &
365	id_var_zw_mask(mid,av) )
366	CALL handle_netcdf_error( 'netcdf', 9998 )
367
368	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
369	'units', 'meters' )
370	CALL handle_netcdf_error( 'netcdf', 9998 )
371
372	!
373	!-- Define x-axis (for scalar position)
374	nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'x', &
375	mask_size(mid,1), id_dim_x_mask(mid,av) )
376	CALL handle_netcdf_error( 'netcdf', 9998 )
377
378	nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'x', NF90_DOUBLE, &
379	id_dim_x_mask(mid,av), id_var_x_mask(mid,av) )
380	CALL handle_netcdf_error( 'netcdf', 9998 )
381
382	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_x_mask(mid,av), &
383	'units', 'meters' )
384	CALL handle_netcdf_error( 'netcdf', 9998 )
385
386	!
387	!-- Define x-axis (for u position)
388	nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'xu', &
389	mask_size(mid,1), id_dim_xu_mask(mid,av) )
390	CALL handle_netcdf_error( 'netcdf', 9998 )
391
392	nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'xu', NF90_DOUBLE, &
393	id_dim_xu_mask(mid,av), &
394	id_var_xu_mask(mid,av) )
395	CALL handle_netcdf_error( 'netcdf', 9998 )
396
397	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_xu_mask(mid,av), &
398	'units', 'meters' )
399	CALL handle_netcdf_error( 'netcdf', 9998 )
400
401	!
402	!-- Define y-axis (for scalar position)
403	nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'y', &
404	mask_size(mid,2), id_dim_y_mask(mid,av) )
405	CALL handle_netcdf_error( 'netcdf', 9998 )
406
407	nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'y', NF90_DOUBLE, &
408	id_dim_y_mask(mid,av), id_var_y_mask(mid,av) )
409	CALL handle_netcdf_error( 'netcdf', 9998 )
410
411	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_y_mask(mid,av), &
412	'units', 'meters' )
413	CALL handle_netcdf_error( 'netcdf', 9998 )
414
415	!
416	!-- Define y-axis (for v position)
417	nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'yv', &
418	mask_size(mid,2), id_dim_yv_mask(mid,av) )
419	CALL handle_netcdf_error( 'netcdf', 9998 )
420
421	nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'yv', NF90_DOUBLE, &
422	id_dim_yv_mask(mid,av), &
423	id_var_yv_mask(mid,av) )
424	CALL handle_netcdf_error( 'netcdf', 9998 )
425
426	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_yv_mask(mid,av), &
427	'units', 'meters' )
428	CALL handle_netcdf_error( 'netcdf', 9998 )
429
430	!
431	!-- In case of non-flat topography define 2d-arrays containing the height
432	!-- informations
433	IF ( TRIM( topography ) /= 'flat' ) THEN
434	!
435	!-- Define zusi = zu(nzb_s_inner)
436	nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zusi', NF90_DOUBLE, &
437	(/ id_dim_x_mask(mid,av), &
438	id_dim_y_mask(mid,av) /), &
439	id_var_zusi_mask(mid,av) )
440	CALL handle_netcdf_error( 'netcdf', 9998 )
441
442	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
443	id_var_zusi_mask(mid,av), &
444	'units', 'meters' )
445	CALL handle_netcdf_error( 'netcdf', 9998 )
446
447	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
448	id_var_zusi_mask(mid,av), &
449	'long_name', 'zu(nzb_s_inner)' )
450	CALL handle_netcdf_error( 'netcdf', 9998 )
451
452	!
453	!-- Define zwwi = zw(nzb_w_inner)
454	nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zwwi', NF90_DOUBLE, &
455	(/ id_dim_x_mask(mid,av), &
456	id_dim_y_mask(mid,av) /), &
457	id_var_zwwi_mask(mid,av) )
458	CALL handle_netcdf_error( 'netcdf', 9998 )
459
460	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
461	id_var_zwwi_mask(mid,av), &
462	'units', 'meters' )
463	CALL handle_netcdf_error( 'netcdf', 9998 )
464
465	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
466	id_var_zwwi_mask(mid,av), &
467	'long_name', 'zw(nzb_w_inner)' )
468	CALL handle_netcdf_error( 'netcdf', 9998 )
469
470	ENDIF
471
472
473	!
474	!-- Define the variables
475	var_list = ';'
476	i = 1
477
478	DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
479
480	!
481	!-- Check for the grid
482	found = .TRUE.
483	SELECT CASE ( domask(mid,av,i) )
484	!
485	!-- Most variables are defined on the scalar grid
486	CASE ( 'e', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c', 'ql_v', &
487	'ql_vp', 'qv', 'rho', 's', 'sa', 'vpt' )
488
489	grid_x = 'x'
490	grid_y = 'y'
491	grid_z = 'zu'
492	!
493	!-- u grid
494	CASE ( 'u' )
495
496	grid_x = 'xu'
497	grid_y = 'y'
498	grid_z = 'zu'
499	!
500	!-- v grid
501	CASE ( 'v' )
502
503	grid_x = 'x'
504	grid_y = 'yv'
505	grid_z = 'zu'
506	!
507	!-- w grid
508	CASE ( 'w' )
509
510	grid_x = 'x'
511	grid_y = 'y'
512	grid_z = 'zw'
513
514	CASE DEFAULT
515	!
516	!-- Check for user-defined quantities
517	CALL user_define_netcdf_grid( domask(mid,av,i), found, &
518	grid_x, grid_y, grid_z )
519
520	END SELECT
521
522	!
523	!-- Select the respective dimension ids
524	IF ( grid_x == 'x' ) THEN
525	id_x = id_dim_x_mask(mid,av)
526	ELSEIF ( grid_x == 'xu' ) THEN
527	id_x = id_dim_xu_mask(mid,av)
528	ENDIF
529
530	IF ( grid_y == 'y' ) THEN
531	id_y = id_dim_y_mask(mid,av)
532	ELSEIF ( grid_y == 'yv' ) THEN
533	id_y = id_dim_yv_mask(mid,av)
534	ENDIF
535
536	IF ( grid_z == 'zu' ) THEN
537	id_z = id_dim_zu_mask(mid,av)
538	ELSEIF ( grid_z == 'zw' ) THEN
539	id_z = id_dim_zw_mask(mid,av)
540	ENDIF
541
542	!
543	!-- Define the grid
544	nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), domask(mid,av,i), &
545	nc_precision(10+mid+av*20), &
546	(/ id_x, id_y, id_z, &
547	id_dim_time_mask(mid,av) /), &
548	id_var_domask(mid,av,i) )
549
550	IF ( .NOT. found ) THEN
551	WRITE ( message_string, * ) 'no grid defined for', &
552	' variable ', TRIM( domask(mid,av,i) )
553	CALL message( 'define_netcdf_header', 'PA9998', 0, 1, 0, 6, 0 )
554	ENDIF
555
556	var_list = TRIM( var_list ) // TRIM( domask(mid,av,i) ) // ';'
557
558	CALL handle_netcdf_error( 'netcdf', 9998 )
559	!
560	!-- Store the 'real' name of the variable (with *, for example)
561	!-- in the long_name attribute. This is evaluated by Ferret,
562	!-- for example.
563	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
564	id_var_domask(mid,av,i), &
565	'long_name', domask(mid,av,i) )
566	CALL handle_netcdf_error( 'netcdf', 9998 )
567	!
568	!-- Define the variable's unit
569	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
570	id_var_domask(mid,av,i), &
571	'units', TRIM( domask_unit(mid,av,i) ) )
572	CALL handle_netcdf_error( 'netcdf', 9998 )
573
574	i = i + 1
575
576	ENDDO
577
578	!
579	!-- No arrays to output
580	IF ( i == 1 ) RETURN
581
582	!
583	!-- Write the list of variables as global attribute (this is used by
584	!-- restart runs and by combine_plot_fields)
585	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
586	'VAR_LIST', var_list )
587	CALL handle_netcdf_error( 'netcdf', 9998 )
588
589	!
590	!-- Leave NetCDF define mode
591	nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
592	CALL handle_netcdf_error( 'netcdf', 9998 )
593
594	!
595	!-- Write data for x (shifted by +dx/2) and xu axis
596	ALLOCATE( netcdf_data(mask_size(mid,1)) )
597
598	netcdf_data = ( mask_i_global(mid,:mask_size(mid,1)) + 0.5 ) * dx
599
600	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_x_mask(mid,av), &
601	netcdf_data, start = (/ 1 /), &
602	count = (/ mask_size(mid,1) /) )
603	CALL handle_netcdf_error( 'netcdf', 9998 )
604
605	netcdf_data = mask_i_global(mid,:mask_size(mid,1)) * dx
606
607	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_xu_mask(mid,av),&
608	netcdf_data, start = (/ 1 /), &
609	count = (/ mask_size(mid,1) /) )
610	CALL handle_netcdf_error( 'netcdf', 9998 )
611
612	DEALLOCATE( netcdf_data )
613
614	!
615	!-- Write data for y (shifted by +dy/2) and yv axis
616	ALLOCATE( netcdf_data(mask_size(mid,2)) )
617
618	netcdf_data = ( mask_j_global(mid,:mask_size(mid,2)) + 0.5 ) * dy
619
620	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_y_mask(mid,av), &
621	netcdf_data, start = (/ 1 /), &
622	count = (/ mask_size(mid,2) /))
623	CALL handle_netcdf_error( 'netcdf', 9998 )
624
625	netcdf_data = mask_j_global(mid,:mask_size(mid,2)) * dy
626
627	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_yv_mask(mid,av), &
628	netcdf_data, start = (/ 1 /), &
629	count = (/ mask_size(mid,2) /))
630	CALL handle_netcdf_error( 'netcdf', 9998 )
631
632	DEALLOCATE( netcdf_data )
633
634	!
635	!-- Write zu and zw data (vertical axes)
636	ALLOCATE( netcdf_data(mask_size(mid,3)) )
637
638	netcdf_data = zu( mask_k_global(mid,:mask_size(mid,3)) )
639
640	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
641	netcdf_data, start = (/ 1 /), &
642	count = (/ mask_size(mid,3) /) )
643	CALL handle_netcdf_error( 'netcdf', 9998 )
644
645	netcdf_data = zw( mask_k_global(mid,:mask_size(mid,3)) )
646
647	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
648	netcdf_data, start = (/ 1 /), &
649	count = (/ mask_size(mid,3) /) )
650	CALL handle_netcdf_error( 'netcdf', 9998 )
651
652	DEALLOCATE( netcdf_data )
653
654	!
655	!-- In case of non-flat topography write height information
656	IF ( TRIM( topography ) /= 'flat' ) THEN
657
658	ALLOCATE( netcdf_data_2d(mask_size(mid,1),mask_size(mid,2)) )
659	netcdf_data_2d = zu_s_inner( mask_i_global(mid,:mask_size(mid,1)),&
660	mask_j_global(mid,:mask_size(mid,2)) )
661
662	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
663	id_var_zusi_mask(mid,av), &
664	netcdf_data_2d, &
665	start = (/ 1, 1 /), &
666	count = (/ mask_size(mid,1), &
667	mask_size(mid,2) /) )
668	CALL handle_netcdf_error( 'netcdf', 9998 )
669
670	netcdf_data_2d = zw_w_inner( mask_i_global(mid,:mask_size(mid,1)),&
671	mask_j_global(mid,:mask_size(mid,2)) )
672
673	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
674	id_var_zwwi_mask(mid,av), &
675	netcdf_data_2d, &
676	start = (/ 1, 1 /), &
677	count = (/ mask_size(mid,1), &
678	mask_size(mid,2) /) )
679	CALL handle_netcdf_error( 'netcdf', 9998 )
680
681	DEALLOCATE( netcdf_data_2d )
682
683	ENDIF
684	!
685	!-- restore original parameter file_id (=formal parameter av) into av
686	av = file_id
687
688
689	CASE ( 'ma_ext' )
690
691	!
692	!-- decompose actual parameter file_id (=formal parameter av) into
693	!-- mid and av
694	file_id = av
695	IF ( file_id <= 140 ) THEN
696	mid = file_id - 120
697	av = 0
698	ELSE
699	mid = file_id - 140
700	av = 1
701	ENDIF
702
703	!
704	!-- Get the list of variables and compare with the actual run.
705	!-- First var_list_old has to be reset, since GET_ATT does not assign
706	!-- trailing blanks.
707	var_list_old = ' '
708	nc_stat = NF90_GET_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'VAR_LIST',&
709	var_list_old )
710	CALL handle_netcdf_error( 'netcdf', 9998 )
711
712	var_list = ';'
713	i = 1
714	DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
715	var_list = TRIM(var_list) // TRIM( domask(mid,av,i) ) // ';'
716	i = i + 1
717	ENDDO
718
719	IF ( av == 0 ) THEN
720	var = '(mask)'
721	ELSE
722	var = '(mask_av)'
723	ENDIF
724
725	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
726	WRITE ( message_string, * ) 'NetCDF file for ', TRIM( var ), &
727	' data for mask', mid, ' from previous run found,', &
728	'&but this file cannot be extended due to variable ', &
729	'mismatch.&New file is created instead.'
730	CALL message( 'define_netcdf_header', 'PA9998', 0, 1, 0, 6, 0 )
731	extend = .FALSE.
732	RETURN
733	ENDIF
734
735	!
736	!-- Get and compare the number of vertical gridpoints
737	nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'zu', &
738	id_var_zu_mask(mid,av) )
739	CALL handle_netcdf_error( 'netcdf', 9998 )
740
741	nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av), &
742	id_var_zu_mask(mid,av), &
743	dimids = id_dim_zu_mask_old )
744	CALL handle_netcdf_error( 'netcdf', 9998 )
745	id_dim_zu_mask(mid,av) = id_dim_zu_mask_old(1)
746
747	nc_stat = NF90_INQUIRE_DIMENSION( id_set_mask(mid,av), &
748	id_dim_zu_mask(mid,av), &
749	len = nz_old )
750	CALL handle_netcdf_error( 'netcdf', 9998 )
751
752	IF ( mask_size(mid,3) /= nz_old ) THEN
753	WRITE ( message_string, * ) 'NetCDF file for ', TRIM( var ), &
754	' data for mask', mid, ' from previous run found,', &
755	'&but this file cannot be extended due to mismatch in ', &
756	' number of&vertical grid points.', &
757	'&New file is created instead.'
758	CALL message( 'define_netcdf_header', 'PA9998', 0, 1, 0, 6, 0 )
759	extend = .FALSE.
760	RETURN
761	ENDIF
762
763	!
764	!-- Get the id of the time coordinate (unlimited coordinate) and its
765	!-- last index on the file. The next time level is plmask..count+1.
766	!-- The current time must be larger than the last output time
767	!-- on the file.
768	nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'time', &
769	id_var_time_mask(mid,av) )
770	CALL handle_netcdf_error( 'netcdf', 9998 )
771
772	nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av), &
773	id_var_time_mask(mid,av), &
774	dimids = id_dim_time_old )
775	CALL handle_netcdf_error( 'netcdf', 9998 )
776	id_dim_time_mask(mid,av) = id_dim_time_old(1)
777
778	nc_stat = NF90_INQUIRE_DIMENSION( id_set_mask(mid,av), &
779	id_dim_time_mask(mid,av), &
780	len = domask_time_count(mid,av) )
781	CALL handle_netcdf_error( 'netcdf', 9998 )
782
783	nc_stat = NF90_GET_VAR( id_set_mask(mid,av), &
784	id_var_time_mask(mid,av), &
785	last_time_coordinate, &
786	start = (/ domask_time_count(mid,av) /), &
787	count = (/ 1 /) )
788	CALL handle_netcdf_error( 'netcdf', 9998 )
789
790	IF ( last_time_coordinate(1) >= simulated_time ) THEN
791	WRITE ( message_string, * ) 'NetCDF file for ', TRIM( var ), &
792	' data for mask', mid, ' from previous run found,', &
793	'&but this file cannot be extended because the current ', &
794	'output time&is less or equal than the last output time ', &
795	'on this file.&New file is created instead.'
796	CALL message( 'define_netcdf_header', 'PA9998', 0, 1, 0, 6, 0 )
797	domask_time_count(mid,av) = 0
798	extend = .FALSE.
799	RETURN
800	ENDIF
801
802	!
803	!-- Dataset seems to be extendable.
804	!-- Now get the variable ids.
805	i = 1
806	DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
807	nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), &
808	TRIM( domask(mid,av,i) ), &
809	id_var_domask(mid,av,i) )
810	CALL handle_netcdf_error( 'netcdf', 9998 )
811	i = i + 1
812	ENDDO
813
814	!
815	!-- Update the title attribute on file
816	!-- In order to avoid 'data mode' errors if updated attributes are larger
817	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
818	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
819	!-- performance loss due to data copying; an alternative strategy would be
820	!-- to ensure equal attribute size in a job chain. Maybe revise later.
821	IF ( av == 0 ) THEN
822	time_average_text = ' '
823	ELSE
824	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
825	averaging_interval
826	ENDIF
827	nc_stat = NF90_REDEF( id_set_mask(mid,av) )
828	CALL handle_netcdf_error( 'netcdf', 9998 )
829	nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'title', &
830	TRIM( run_description_header ) // &
831	TRIM( time_average_text ) )
832	CALL handle_netcdf_error( 'netcdf', 9998 )
833	nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
834	CALL handle_netcdf_error( 'netcdf', 9998 )
835	WRITE ( message_string, * ) 'NetCDF file for ', TRIM( var ), &
836	' data for mask', mid, ' from previous run found.', &
837	'&This file will be extended.'
838	CALL message( 'define_netcdf_header', 'PA9998', 0, 0, 0, 6, 0 )
839	!
840	!-- restore original parameter file_id (=formal parameter av) into av
841	av = file_id
842
843
844	CASE ( '3d_new' )
845
846	!
847	!-- Define some global attributes of the dataset
848	nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'Conventions', &
849	'COARDS' )
850	CALL handle_netcdf_error( 'netcdf', 62 )
851
852	IF ( av == 0 ) THEN
853	time_average_text = ' '
854	ELSE
855	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
856	averaging_interval
857	ENDIF
858	nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'title', &
859	TRIM( run_description_header ) // &
860	TRIM( time_average_text ) )
861	CALL handle_netcdf_error( 'netcdf', 63 )
862	IF ( av == 1 ) THEN
863	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval
864	nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'time_avg', &
865	TRIM( time_average_text ) )
866	CALL handle_netcdf_error( 'netcdf', 63 )
867	ENDIF
868
869	!
870	!-- Define time coordinate for volume data (unlimited dimension)
871	nc_stat = NF90_DEF_DIM( id_set_3d(av), 'time', NF90_UNLIMITED, &
872	id_dim_time_3d(av) )
873	CALL handle_netcdf_error( 'netcdf', 64 )
874
875	nc_stat = NF90_DEF_VAR( id_set_3d(av), 'time', NF90_DOUBLE, &
876	id_dim_time_3d(av), id_var_time_3d(av) )
877	CALL handle_netcdf_error( 'netcdf', 65 )
878
879	nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_time_3d(av), 'units', &
880	'seconds')
881	CALL handle_netcdf_error( 'netcdf', 66 )
882
883	!
884	!-- Define spatial dimensions and coordinates:
885	!-- Define vertical coordinate grid (zu grid)
886	nc_stat = NF90_DEF_DIM( id_set_3d(av), 'zu_3d', nz_do3d-nzb+1, &
887	id_dim_zu_3d(av) )
888	CALL handle_netcdf_error( 'netcdf', 67 )
889
890	nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zu_3d', NF90_DOUBLE, &
891	id_dim_zu_3d(av), id_var_zu_3d(av) )
892	CALL handle_netcdf_error( 'netcdf', 68 )
893
894	nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zu_3d(av), 'units', &
895	'meters' )
896	CALL handle_netcdf_error( 'netcdf', 69 )
897
898	!
899	!-- Define vertical coordinate grid (zw grid)
900	nc_stat = NF90_DEF_DIM( id_set_3d(av), 'zw_3d', nz_do3d-nzb+1, &
901	id_dim_zw_3d(av) )
902	CALL handle_netcdf_error( 'netcdf', 70 )
903
904	nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zw_3d', NF90_DOUBLE, &
905	id_dim_zw_3d(av), id_var_zw_3d(av) )
906	CALL handle_netcdf_error( 'netcdf', 71 )
907
908	nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zw_3d(av), 'units', &
909	'meters' )
910	CALL handle_netcdf_error( 'netcdf', 72 )
911
912	!
913	!-- Define x-axis (for scalar position)
914	nc_stat = NF90_DEF_DIM( id_set_3d(av), 'x', nx+2, id_dim_x_3d(av) )
915	CALL handle_netcdf_error( 'netcdf', 73 )
916
917	nc_stat = NF90_DEF_VAR( id_set_3d(av), 'x', NF90_DOUBLE, &
918	id_dim_x_3d(av), id_var_x_3d(av) )
919	CALL handle_netcdf_error( 'netcdf', 74 )
920
921	nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_x_3d(av), 'units', &
922	'meters' )
923	CALL handle_netcdf_error( 'netcdf', 75 )
924
925	!
926	!-- Define x-axis (for u position)
927	nc_stat = NF90_DEF_DIM( id_set_3d(av), 'xu', nx+2, id_dim_xu_3d(av) )
928	CALL handle_netcdf_error( 'netcdf', 358 )
929
930	nc_stat = NF90_DEF_VAR( id_set_3d(av), 'xu', NF90_DOUBLE, &
931	id_dim_xu_3d(av), id_var_xu_3d(av) )
932	CALL handle_netcdf_error( 'netcdf', 359 )
933
934	nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_xu_3d(av), 'units', &
935	'meters' )
936	CALL handle_netcdf_error( 'netcdf', 360 )
937
938	!
939	!-- Define y-axis (for scalar position)
940	nc_stat = NF90_DEF_DIM( id_set_3d(av), 'y', ny+2, id_dim_y_3d(av) )
941	CALL handle_netcdf_error( 'netcdf', 76 )
942
943	nc_stat = NF90_DEF_VAR( id_set_3d(av), 'y', NF90_DOUBLE, &
944	id_dim_y_3d(av), id_var_y_3d(av) )
945	CALL handle_netcdf_error( 'netcdf', 77 )
946
947	nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_y_3d(av), 'units', &
948	'meters' )
949	CALL handle_netcdf_error( 'netcdf', 78 )
950
951	!
952	!-- Define y-axis (for v position)
953	nc_stat = NF90_DEF_DIM( id_set_3d(av), 'yv', ny+2, id_dim_yv_3d(av) )
954	CALL handle_netcdf_error( 'netcdf', 361 )
955
956	nc_stat = NF90_DEF_VAR( id_set_3d(av), 'yv', NF90_DOUBLE, &
957	id_dim_yv_3d(av), id_var_yv_3d(av) )
958	CALL handle_netcdf_error( 'netcdf', 362 )
959
960	nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_yv_3d(av), 'units', &
961	'meters' )
962	CALL handle_netcdf_error( 'netcdf', 363 )
963
964	!
965	!-- In case of non-flat topography define 2d-arrays containing the height
966	!-- informations
967	IF ( TRIM( topography ) /= 'flat' ) THEN
968	!
969	!-- Define zusi = zu(nzb_s_inner)
970	nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zusi', NF90_DOUBLE, &
971	(/ id_dim_x_3d(av), id_dim_y_3d(av) /), &
972	id_var_zusi_3d(av) )
973	CALL handle_netcdf_error( 'netcdf', 413 )
974
975	nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zusi_3d(av), &
976	'units', 'meters' )
977	CALL handle_netcdf_error( 'netcdf', 414 )
978
979	nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zusi_3d(av), &
980	'long_name', 'zu(nzb_s_inner)' )
981	CALL handle_netcdf_error( 'netcdf', 415 )
982
983	!
984	!-- Define zwwi = zw(nzb_w_inner)
985	nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zwwi', NF90_DOUBLE, &
986	(/ id_dim_x_3d(av), id_dim_y_3d(av) /), &
987	id_var_zwwi_3d(av) )
988	CALL handle_netcdf_error( 'netcdf', 416 )
989
990	nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zwwi_3d(av), &
991	'units', 'meters' )
992	CALL handle_netcdf_error( 'netcdf', 417 )
993
994	nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zwwi_3d(av), &
995	'long_name', 'zw(nzb_w_inner)' )
996	CALL handle_netcdf_error( 'netcdf', 418 )
997
998	ENDIF
999
1000
1001	!
1002	!-- Define the variables
1003	var_list = ';'
1004	i = 1
1005
1006	DO WHILE ( do3d(av,i)(1:1) /= ' ' )
1007
1008	!
1009	!-- Check for the grid
1010	found = .TRUE.
1011	SELECT CASE ( do3d(av,i) )
1012	!
1013	!-- Most variables are defined on the scalar grid
1014	CASE ( 'e', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c', 'ql_v', &
1015	'ql_vp', 'qv', 'rho', 's', 'sa', 'vpt' )
1016
1017	grid_x = 'x'
1018	grid_y = 'y'
1019	grid_z = 'zu'
1020	!
1021	!-- u grid
1022	CASE ( 'u' )
1023
1024	grid_x = 'xu'
1025	grid_y = 'y'
1026	grid_z = 'zu'
1027	!
1028	!-- v grid
1029	CASE ( 'v' )
1030
1031	grid_x = 'x'
1032	grid_y = 'yv'
1033	grid_z = 'zu'
1034	!
1035	!-- w grid
1036	CASE ( 'w' )
1037
1038	grid_x = 'x'
1039	grid_y = 'y'
1040	grid_z = 'zw'
1041
1042	CASE DEFAULT
1043	!
1044	!-- Check for user-defined quantities
1045	CALL user_define_netcdf_grid( do3d(av,i), found, grid_x, &
1046	grid_y, grid_z )
1047
1048	END SELECT
1049
1050	!
1051	!-- Select the respective dimension ids
1052	IF ( grid_x == 'x' ) THEN
1053	id_x = id_dim_x_3d(av)
1054	ELSEIF ( grid_x == 'xu' ) THEN
1055	id_x = id_dim_xu_3d(av)
1056	ENDIF
1057
1058	IF ( grid_y == 'y' ) THEN
1059	id_y = id_dim_y_3d(av)
1060	ELSEIF ( grid_y == 'yv' ) THEN
1061	id_y = id_dim_yv_3d(av)
1062	ENDIF
1063
1064	IF ( grid_z == 'zu' ) THEN
1065	id_z = id_dim_zu_3d(av)
1066	ELSEIF ( grid_z == 'zw' ) THEN
1067	id_z = id_dim_zw_3d(av)
1068	ENDIF
1069
1070	!
1071	!-- Define the grid
1072	nc_stat = NF90_DEF_VAR( id_set_3d(av), do3d(av,i), &
1073	nc_precision(4), &
1074	(/ id_x, id_y, id_z, id_dim_time_3d(av) /), &
1075	id_var_do3d(av,i) )
1076
1077	IF ( .NOT. found ) THEN
1078	WRITE ( message_string, * ) 'no grid defined for', &
1079	' variable ', TRIM( do3d(av,i) )
1080	CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, 6, 0 )
1081	ENDIF
1082
1083	var_list = TRIM( var_list ) // TRIM( do3d(av,i) ) // ';'
1084
1085	CALL handle_netcdf_error( 'netcdf', 79 )
1086	!
1087	!-- Store the 'real' name of the variable (with *, for example)
1088	!-- in the long_name attribute. This is evaluated by Ferret,
1089	!-- for example.
1090	nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_do3d(av,i), &
1091	'long_name', do3d(av,i) )
1092	CALL handle_netcdf_error( 'netcdf', 80 )
1093	!
1094	!-- Define the variable's unit
1095	nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_do3d(av,i), &
1096	'units', TRIM( do3d_unit(av,i) ) )
1097	CALL handle_netcdf_error( 'netcdf', 357 )
1098	#if defined( __netcdf4 )
1099	!
1100	!-- Set collective io operations for parallel io
1101	IF ( netcdf_data_format > 2 ) THEN
1102	nc_stat = NF90_VAR_PAR_ACCESS( id_set_3d(av), &
1103	id_var_do3d(av,i), &
1104	NF90_COLLECTIVE )
1105	CALL handle_netcdf_error( 'netcdf', 445 )
1106	ENDIF
1107	#endif
1108	i = i + 1
1109
1110	ENDDO
1111
1112	!
1113	!-- No arrays to output
1114	IF ( i == 1 ) RETURN
1115
1116	!
1117	!-- Write the list of variables as global attribute (this is used by
1118	!-- restart runs and by combine_plot_fields)
1119	nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
1120	var_list )
1121	CALL handle_netcdf_error( 'netcdf', 81 )
1122
1123	!
1124	!-- Leave NetCDF define mode
1125	nc_stat = NF90_ENDDEF( id_set_3d(av) )
1126	CALL handle_netcdf_error( 'netcdf', 82 )
1127
1128	!
1129	!-- Write data for x (shifted by +dx/2) and xu axis
1130	ALLOCATE( netcdf_data(0:nx+1) )
1131
1132	DO i = 0, nx+1
1133	netcdf_data(i) = ( i + 0.5 ) * dx
1134	ENDDO
1135
1136	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_x_3d(av), netcdf_data, &
1137	start = (/ 1 /), count = (/ nx+2 /) )
1138	CALL handle_netcdf_error( 'netcdf', 83 )
1139
1140	DO i = 0, nx+1
1141	netcdf_data(i) = i * dx
1142	ENDDO
1143
1144	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_xu_3d(av), &
1145	netcdf_data, start = (/ 1 /), &
1146	count = (/ nx+2 /) )
1147	CALL handle_netcdf_error( 'netcdf', 385 )
1148
1149	DEALLOCATE( netcdf_data )
1150
1151	!
1152	!-- Write data for y (shifted by +dy/2) and yv axis
1153	ALLOCATE( netcdf_data(0:ny+1) )
1154
1155	DO i = 0, ny+1
1156	netcdf_data(i) = ( i + 0.5 ) * dy
1157	ENDDO
1158
1159	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_y_3d(av), netcdf_data, &
1160	start = (/ 1 /), count = (/ ny+2 /))
1161	CALL handle_netcdf_error( 'netcdf', 84 )
1162
1163	DO i = 0, ny+1
1164	netcdf_data(i) = i * dy
1165	ENDDO
1166
1167	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_yv_3d(av), &
1168	netcdf_data, start = (/ 1 /), &
1169	count = (/ ny+2 /))
1170	CALL handle_netcdf_error( 'netcdf', 387 )
1171
1172	DEALLOCATE( netcdf_data )
1173
1174	!
1175	!-- Write zu and zw data (vertical axes)
1176	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zu_3d(av), &
1177	zu(nzb:nz_do3d), start = (/ 1 /), &
1178	count = (/ nz_do3d-nzb+1 /) )
1179	CALL handle_netcdf_error( 'netcdf', 85 )
1180
1181
1182	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zw_3d(av), &
1183	zw(nzb:nz_do3d), start = (/ 1 /), &
1184	count = (/ nz_do3d-nzb+1 /) )
1185	CALL handle_netcdf_error( 'netcdf', 86 )
1186
1187
1188	!
1189	!-- In case of non-flat topography write height information
1190	IF ( TRIM( topography ) /= 'flat' ) THEN
1191
1192	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av), &
1193	zu_s_inner(0:nx+1,0:ny+1), &
1194	start = (/ 1, 1 /), &
1195	count = (/ nx+2, ny+2 /) )
1196	CALL handle_netcdf_error( 'netcdf', 419 )
1197
1198	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av), &
1199	zw_w_inner(0:nx+1,0:ny+1), &
1200	start = (/ 1, 1 /), &
1201	count = (/ nx+2, ny+2 /) )
1202	CALL handle_netcdf_error( 'netcdf', 420 )
1203
1204	ENDIF
1205
1206	CASE ( '3d_ext' )
1207
1208	!
1209	!-- Get the list of variables and compare with the actual run.
1210	!-- First var_list_old has to be reset, since GET_ATT does not assign
1211	!-- trailing blanks.
1212	var_list_old = ' '
1213	nc_stat = NF90_GET_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
1214	var_list_old )
1215	CALL handle_netcdf_error( 'netcdf', 87 )
1216
1217	var_list = ';'
1218	i = 1
1219	DO WHILE ( do3d(av,i)(1:1) /= ' ' )
1220	var_list = TRIM(var_list) // TRIM( do3d(av,i) ) // ';'
1221	i = i + 1
1222	ENDDO
1223
1224	IF ( av == 0 ) THEN
1225	var = '(3d)'
1226	ELSE
1227	var = '(3d_av)'
1228	ENDIF
1229
1230	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
1231	message_string = 'NetCDF file for volume data ' // &
1232	TRIM( var ) // ' from previous run found,' // &
1233	'&but this file cannot be extended due to' // &
1234	' variable mismatch.' // &
1235	'&New file is created instead.'
1236	CALL message( 'define_netcdf_header', 'PA0245', 0, 1, 0, 6, 0 )
1237	extend = .FALSE.
1238	RETURN
1239	ENDIF
1240
1241	!
1242	!-- Get and compare the number of vertical gridpoints
1243	nc_stat = NF90_INQ_VARID( id_set_3d(av), 'zu_3d', id_var_zu_3d(av) )
1244	CALL handle_netcdf_error( 'netcdf', 88 )
1245
1246	nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_zu_3d(av), &
1247	dimids = id_dim_zu_3d_old )
1248	CALL handle_netcdf_error( 'netcdf', 89 )
1249	id_dim_zu_3d(av) = id_dim_zu_3d_old(1)
1250
1251	nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_zu_3d(av), &
1252	len = nz_old )
1253	CALL handle_netcdf_error( 'netcdf', 90 )
1254
1255	IF ( nz_do3d-nzb+1 /= nz_old ) THEN
1256	message_string = 'NetCDF file for volume data ' // &
1257	TRIM( var ) // ' from previous run found,' // &
1258	'&but this file cannot be extended due to' // &
1259	' mismatch in number of' // &
1260	'&vertical grid points (nz_do3d).' // &
1261	'&New file is created instead.'
1262	CALL message( 'define_netcdf_header', 'PA0246', 0, 1, 0, 6, 0 )
1263	extend = .FALSE.
1264	RETURN
1265	ENDIF
1266
1267	!
1268	!-- Get the id of the time coordinate (unlimited coordinate) and its
1269	!-- last index on the file. The next time level is pl3d..count+1.
1270	!-- The current time must be larger than the last output time
1271	!-- on the file.
1272	nc_stat = NF90_INQ_VARID( id_set_3d(av), 'time', id_var_time_3d(av) )
1273	CALL handle_netcdf_error( 'netcdf', 91 )
1274
1275	nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_time_3d(av), &
1276	dimids = id_dim_time_old )
1277	CALL handle_netcdf_error( 'netcdf', 92 )
1278
1279	id_dim_time_3d(av) = id_dim_time_old(1)
1280
1281	nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_time_3d(av), &
1282	len = do3d_time_count(av) )
1283	CALL handle_netcdf_error( 'netcdf', 93 )
1284
1285	nc_stat = NF90_GET_VAR( id_set_3d(av), id_var_time_3d(av), &
1286	last_time_coordinate, &
1287	start = (/ do3d_time_count(av) /), &
1288	count = (/ 1 /) )
1289	CALL handle_netcdf_error( 'netcdf', 94 )
1290
1291	IF ( last_time_coordinate(1) >= simulated_time ) THEN
1292	message_string = 'NetCDF file for volume data ' // &
1293	TRIM( var ) // ' from previous run found,' // &
1294	'&but this file cannot be extended becaus' // &
1295	'e the current output time' // &
1296	'&is less or equal than the last output t' // &
1297	'ime on this file.' // &
1298	'&New file is created instead.'
1299	CALL message( 'define_netcdf_header', 'PA0247', 0, 1, 0, 6, 0 )
1300	do3d_time_count(av) = 0
1301	extend = .FALSE.
1302	RETURN
1303	ENDIF
1304
1305	!
1306	!-- Dataset seems to be extendable.
1307	!-- Now get the variable ids.
1308	i = 1
1309	DO WHILE ( do3d(av,i)(1:1) /= ' ' )
1310	nc_stat = NF90_INQ_VARID( id_set_3d(av), TRIM( do3d(av,i) ), &
1311	id_var_do3d(av,i) )
1312	CALL handle_netcdf_error( 'netcdf', 95 )
1313	#if defined( __netcdf4 )
1314	!
1315	!-- Set collective io operations for parallel io
1316	IF ( netcdf_data_format > 2 ) THEN
1317	nc_stat = NF90_VAR_PAR_ACCESS( id_set_3d(av), &
1318	id_var_do3d(av,i), &
1319	NF90_COLLECTIVE )
1320	CALL handle_netcdf_error( 'netcdf', 453 )
1321	ENDIF
1322	#endif
1323	i = i + 1
1324	ENDDO
1325
1326	!
1327	!-- Update the title attribute on file
1328	!-- In order to avoid 'data mode' errors if updated attributes are larger
1329	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
1330	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
1331	!-- performance loss due to data copying; an alternative strategy would be
1332	!-- to ensure equal attribute size. Maybe revise later.
1333	IF ( av == 0 ) THEN
1334	time_average_text = ' '
1335	ELSE
1336	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
1337	averaging_interval
1338	ENDIF
1339	nc_stat = NF90_REDEF( id_set_3d(av) )
1340	CALL handle_netcdf_error( 'netcdf', 429 )
1341	nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'title', &
1342	TRIM( run_description_header ) // &
1343	TRIM( time_average_text ) )
1344	CALL handle_netcdf_error( 'netcdf', 96 )
1345	nc_stat = NF90_ENDDEF( id_set_3d(av) )
1346	CALL handle_netcdf_error( 'netcdf', 430 )
1347	message_string = 'NetCDF file for volume data ' // &
1348	TRIM( var ) // ' from previous run found.' // &
1349	'&This file will be extended.'
1350	CALL message( 'define_netcdf_header', 'PA0248', 0, 0, 0, 6, 0 )
1351
1352	CASE ( 'xy_new' )
1353
1354	!
1355	!-- Define some global attributes of the dataset
1356	nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'Conventions', &
1357	'COARDS' )
1358	CALL handle_netcdf_error( 'netcdf', 97 )
1359
1360	IF ( av == 0 ) THEN
1361	time_average_text = ' '
1362	ELSE
1363	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
1364	averaging_interval
1365	ENDIF
1366	nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'title', &
1367	TRIM( run_description_header ) // &
1368	TRIM( time_average_text ) )
1369	CALL handle_netcdf_error( 'netcdf', 98 )
1370	IF ( av == 1 ) THEN
1371	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval
1372	nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'time_avg', &
1373	TRIM( time_average_text ) )
1374	CALL handle_netcdf_error( 'netcdf', 98 )
1375	ENDIF
1376
1377	!
1378	!-- Define time coordinate for xy sections (unlimited dimension)
1379	nc_stat = NF90_DEF_DIM( id_set_xy(av), 'time', NF90_UNLIMITED, &
1380	id_dim_time_xy(av) )
1381	CALL handle_netcdf_error( 'netcdf', 99 )
1382
1383	nc_stat = NF90_DEF_VAR( id_set_xy(av), 'time', NF90_DOUBLE, &
1384	id_dim_time_xy(av), id_var_time_xy(av) )
1385	CALL handle_netcdf_error( 'netcdf', 100 )
1386
1387	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_time_xy(av), 'units', &
1388	'seconds')
1389	CALL handle_netcdf_error( 'netcdf', 101 )
1390
1391	!
1392	!-- Define the spatial dimensions and coordinates for xy-sections.
1393	!-- First, determine the number of horizontal sections to be written.
1394	IF ( section(1,1) == -9999 ) THEN
1395	RETURN
1396	ELSE
1397	ns = 1
1398	DO WHILE ( section(ns,1) /= -9999 .AND. ns <= 100 )
1399	ns = ns + 1
1400	ENDDO
1401	ns = ns - 1
1402	ENDIF
1403
1404	!
1405	!-- Define vertical coordinate grid (zu grid)
1406	nc_stat = NF90_DEF_DIM( id_set_xy(av), 'zu_xy', ns, id_dim_zu_xy(av) )
1407	CALL handle_netcdf_error( 'netcdf', 102 )
1408
1409	nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zu_xy', NF90_DOUBLE, &
1410	id_dim_zu_xy(av), id_var_zu_xy(av) )
1411	CALL handle_netcdf_error( 'netcdf', 103 )
1412
1413	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zu_xy(av), 'units', &
1414	'meters' )
1415	CALL handle_netcdf_error( 'netcdf', 104 )
1416
1417	!
1418	!-- Define vertical coordinate grid (zw grid)
1419	nc_stat = NF90_DEF_DIM( id_set_xy(av), 'zw_xy', ns, id_dim_zw_xy(av) )
1420	CALL handle_netcdf_error( 'netcdf', 105 )
1421
1422	nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zw_xy', NF90_DOUBLE, &
1423	id_dim_zw_xy(av), id_var_zw_xy(av) )
1424	CALL handle_netcdf_error( 'netcdf', 106 )
1425
1426	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zw_xy(av), 'units', &
1427	'meters' )
1428	CALL handle_netcdf_error( 'netcdf', 107 )
1429
1430	!
1431	!-- Define a pseudo vertical coordinate grid for the surface variables
1432	!-- u* and t* to store their height level
1433	nc_stat = NF90_DEF_DIM( id_set_xy(av), 'zu1_xy', 1, &
1434	id_dim_zu1_xy(av) )
1435	CALL handle_netcdf_error( 'netcdf', 108 )
1436
1437	nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zu1_xy', NF90_DOUBLE, &
1438	id_dim_zu1_xy(av), id_var_zu1_xy(av) )
1439	CALL handle_netcdf_error( 'netcdf', 109 )
1440
1441	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zu1_xy(av), 'units', &
1442	'meters' )
1443	CALL handle_netcdf_error( 'netcdf', 110 )
1444
1445	!
1446	!-- Define a variable to store the layer indices of the horizontal cross
1447	!-- sections, too
1448	nc_stat = NF90_DEF_VAR( id_set_xy(av), 'ind_z_xy', NF90_DOUBLE, &
1449	id_dim_zu_xy(av), id_var_ind_z_xy(av) )
1450	CALL handle_netcdf_error( 'netcdf', 111 )
1451
1452	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_ind_z_xy(av), 'units', &
1453	'gridpoints')
1454	CALL handle_netcdf_error( 'netcdf', 112 )
1455
1456	!
1457	!-- Define x-axis (for scalar position)
1458	nc_stat = NF90_DEF_DIM( id_set_xy(av), 'x', nx+2, id_dim_x_xy(av) )
1459	CALL handle_netcdf_error( 'netcdf', 113 )
1460
1461	nc_stat = NF90_DEF_VAR( id_set_xy(av), 'x', NF90_DOUBLE, &
1462	id_dim_x_xy(av), id_var_x_xy(av) )
1463	CALL handle_netcdf_error( 'netcdf', 114 )
1464
1465
1466	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_x_xy(av), 'units', &
1467	'meters' )
1468	CALL handle_netcdf_error( 'netcdf', 115 )
1469
1470	!
1471	!-- Define x-axis (for u position)
1472	nc_stat = NF90_DEF_DIM( id_set_xy(av), 'xu', nx+2, id_dim_xu_xy(av) )
1473	CALL handle_netcdf_error( 'netcdf', 388 )
1474
1475	nc_stat = NF90_DEF_VAR( id_set_xy(av), 'xu', NF90_DOUBLE, &
1476	id_dim_xu_xy(av), id_var_xu_xy(av) )
1477	CALL handle_netcdf_error( 'netcdf', 389 )
1478
1479	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_xu_xy(av), 'units', &
1480	'meters' )
1481	CALL handle_netcdf_error( 'netcdf', 390 )
1482
1483	!
1484	!-- Define y-axis (for scalar position)
1485	nc_stat = NF90_DEF_DIM( id_set_xy(av), 'y', ny+2, id_dim_y_xy(av) )
1486	CALL handle_netcdf_error( 'netcdf', 116 )
1487
1488	nc_stat = NF90_DEF_VAR( id_set_xy(av), 'y', NF90_DOUBLE, &
1489	id_dim_y_xy(av), id_var_y_xy(av) )
1490	CALL handle_netcdf_error( 'netcdf', 117 )
1491
1492	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_y_xy(av), 'units', &
1493	'meters' )
1494	CALL handle_netcdf_error( 'netcdf', 118 )
1495
1496	!
1497	!-- Define y-axis (for scalar position)
1498	nc_stat = NF90_DEF_DIM( id_set_xy(av), 'yv', ny+2, id_dim_yv_xy(av) )
1499	CALL handle_netcdf_error( 'netcdf', 364 )
1500
1501	nc_stat = NF90_DEF_VAR( id_set_xy(av), 'yv', NF90_DOUBLE, &
1502	id_dim_yv_xy(av), id_var_yv_xy(av) )
1503	CALL handle_netcdf_error( 'netcdf', 365 )
1504
1505	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_yv_xy(av), 'units', &
1506	'meters' )
1507	CALL handle_netcdf_error( 'netcdf', 366 )
1508
1509	!
1510	!-- In case of non-flat topography define 2d-arrays containing the height
1511	!-- informations
1512	IF ( TRIM( topography ) /= 'flat' ) THEN
1513	!
1514	!-- Define zusi = zu(nzb_s_inner)
1515	nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zusi', NF90_DOUBLE, &
1516	(/ id_dim_x_xy(av), id_dim_y_xy(av) /), &
1517	id_var_zusi_xy(av) )
1518	CALL handle_netcdf_error( 'netcdf', 421 )
1519
1520	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zusi_xy(av), &
1521	'units', 'meters' )
1522	CALL handle_netcdf_error( 'netcdf', 422 )
1523
1524	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zusi_xy(av), &
1525	'long_name', 'zu(nzb_s_inner)' )
1526	CALL handle_netcdf_error( 'netcdf', 423 )
1527
1528	!
1529	!-- Define zwwi = zw(nzb_w_inner)
1530	nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zwwi', NF90_DOUBLE, &
1531	(/ id_dim_x_xy(av), id_dim_y_xy(av) /), &
1532	id_var_zwwi_xy(av) )
1533	CALL handle_netcdf_error( 'netcdf', 424 )
1534
1535	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zwwi_xy(av), &
1536	'units', 'meters' )
1537	CALL handle_netcdf_error( 'netcdf', 425 )
1538
1539	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zwwi_xy(av), &
1540	'long_name', 'zw(nzb_w_inner)' )
1541	CALL handle_netcdf_error( 'netcdf', 426 )
1542
1543	ENDIF
1544
1545
1546	!
1547	!-- Define the variables
1548	var_list = ';'
1549	i = 1
1550
1551	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
1552
1553	IF ( INDEX( do2d(av,i), 'xy' ) /= 0 ) THEN
1554	!
1555	!-- If there is a star in the variable name (u* or t*), it is a
1556	!-- surface variable. Define it with id_dim_zu1_xy.
1557	IF ( INDEX( do2d(av,i), '*' ) /= 0 ) THEN
1558	!
1559	!-- First, remove those characters not allowed by NetCDF
1560	netcdf_var_name = do2d(av,i)
1561	CALL clean_netcdf_varname( netcdf_var_name )
1562
1563	nc_stat = NF90_DEF_VAR( id_set_xy(av), netcdf_var_name, &
1564	nc_precision(1), &
1565	(/ id_dim_x_xy(av), id_dim_y_xy(av),&
1566	id_dim_zu1_xy(av), &
1567	id_dim_time_xy(av) /), &
1568	id_var_do2d(av,i) )
1569
1570	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
1571
1572	ELSE
1573
1574	!
1575	!-- Check for the grid
1576	found = .TRUE.
1577	SELECT CASE ( do2d(av,i) )
1578	!
1579	!-- Most variables are defined on the zu grid
1580	CASE ( 'e_xy', 'p_xy', 'pc_xy', 'pr_xy', 'pt_xy', 'q_xy',&
1581	'ql_xy', 'ql_c_xy', 'ql_v_xy', 'ql_vp_xy', &
1582	'qv_xy', 'rho_xy', 's_xy', 'sa_xy', 'vpt_xy' )
1583
1584	grid_x = 'x'
1585	grid_y = 'y'
1586	grid_z = 'zu'
1587	!
1588	!-- u grid
1589	CASE ( 'u_xy' )
1590
1591	grid_x = 'xu'
1592	grid_y = 'y'
1593	grid_z = 'zu'
1594	!
1595	!-- v grid
1596	CASE ( 'v_xy' )
1597
1598	grid_x = 'x'
1599	grid_y = 'yv'
1600	grid_z = 'zu'
1601	!
1602	!-- w grid
1603	CASE ( 'w_xy' )
1604
1605	grid_x = 'x'
1606	grid_y = 'y'
1607	grid_z = 'zw'
1608
1609	CASE DEFAULT
1610	!
1611	!-- Check for user-defined quantities
1612	CALL user_define_netcdf_grid( do2d(av,i), found, &
1613	grid_x, grid_y, grid_z )
1614
1615	END SELECT
1616
1617	!
1618	!-- Select the respective dimension ids
1619	IF ( grid_x == 'x' ) THEN
1620	id_x = id_dim_x_xy(av)
1621	ELSEIF ( grid_x == 'xu' ) THEN
1622	id_x = id_dim_xu_xy(av)
1623	ENDIF
1624
1625	IF ( grid_y == 'y' ) THEN
1626	id_y = id_dim_y_xy(av)
1627	ELSEIF ( grid_y == 'yv' ) THEN
1628	id_y = id_dim_yv_xy(av)
1629	ENDIF
1630
1631	IF ( grid_z == 'zu' ) THEN
1632	id_z = id_dim_zu_xy(av)
1633	ELSEIF ( grid_z == 'zw' ) THEN
1634	id_z = id_dim_zw_xy(av)
1635	ENDIF
1636
1637	!
1638	!-- Define the grid
1639	nc_stat = NF90_DEF_VAR( id_set_xy(av), do2d(av,i), &
1640	nc_precision(1), &
1641	(/ id_x, id_y, id_z, id_dim_time_xy(av) /), &
1642	id_var_do2d(av,i) )
1643
1644	IF ( .NOT. found ) THEN
1645	WRITE ( message_string, * ) 'no grid defined for', &
1646	' variable ', TRIM( do2d(av,i) )
1647	CALL message( 'define_netcdf_header', 'PA0244',&
1648	0, 1, 0, 6, 0 )
1649	ENDIF
1650
1651	var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
1652
1653	ENDIF
1654
1655	CALL handle_netcdf_error( 'netcdf', 119 )
1656	!
1657	!-- Store the 'real' name of the variable (with *, for example)
1658	!-- in the long_name attribute. This is evaluated by Ferret,
1659	!-- for example.
1660	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_do2d(av,i), &
1661	'long_name', do2d(av,i) )
1662	CALL handle_netcdf_error( 'netcdf', 120 )
1663	!
1664	!-- Define the variable's unit
1665	nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_do2d(av,i), &
1666	'units', TRIM( do2d_unit(av,i) ) )
1667	CALL handle_netcdf_error( 'netcdf', 354 )
1668	#if defined( __netcdf4 )
1669	!
1670	!-- Set collective io operations for parallel io
1671	IF ( netcdf_data_format > 2 ) THEN
1672	nc_stat = NF90_VAR_PAR_ACCESS( id_set_xy(av), &
1673	id_var_do2d(av,i), &
1674	NF90_COLLECTIVE )
1675	CALL handle_netcdf_error( 'netcdf', 448 )
1676	ENDIF
1677	#endif
1678	ENDIF
1679
1680	i = i + 1
1681
1682	ENDDO
1683
1684	!
1685	!-- No arrays to output. Close the netcdf file and return.
1686	IF ( i == 1 ) RETURN
1687
1688	!
1689	!-- Write the list of variables as global attribute (this is used by
1690	!-- restart runs and by combine_plot_fields)
1691	nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'VAR_LIST', &
1692	var_list )
1693	CALL handle_netcdf_error( 'netcdf', 121 )
1694
1695	!
1696	!-- Leave NetCDF define mode
1697	nc_stat = NF90_ENDDEF( id_set_xy(av) )
1698	CALL handle_netcdf_error( 'netcdf', 122 )
1699
1700	!
1701	!-- Write axis data: z_xy, x, y
1702	ALLOCATE( netcdf_data(1:ns) )
1703
1704	!
1705	!-- Write zu data
1706	DO i = 1, ns
1707	IF( section(i,1) == -1 ) THEN
1708	netcdf_data(i) = -1.0 ! section averaged along z
1709	ELSE
1710	netcdf_data(i) = zu( section(i,1) )
1711	ENDIF
1712	ENDDO
1713	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu_xy(av), &
1714	netcdf_data, start = (/ 1 /), &
1715	count = (/ ns /) )
1716	CALL handle_netcdf_error( 'netcdf', 123 )
1717
1718	!
1719	!-- Write zw data
1720	DO i = 1, ns
1721	IF( section(i,1) == -1 ) THEN
1722	netcdf_data(i) = -1.0 ! section averaged along z
1723	ELSE
1724	netcdf_data(i) = zw( section(i,1) )
1725	ENDIF
1726	ENDDO
1727	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zw_xy(av), &
1728	netcdf_data, start = (/ 1 /), &
1729	count = (/ ns /) )
1730	CALL handle_netcdf_error( 'netcdf', 124 )
1731
1732	!
1733	!-- Write gridpoint number data
1734	netcdf_data(1:ns) = section(1:ns,1)
1735	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_ind_z_xy(av), &
1736	netcdf_data, start = (/ 1 /), &
1737	count = (/ ns /) )
1738	CALL handle_netcdf_error( 'netcdf', 125 )
1739
1740	DEALLOCATE( netcdf_data )
1741
1742	!
1743	!-- Write the cross section height u, t
1744	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu1_xy(av), &
1745	(/ zu(nzb+1) /), start = (/ 1 /), &
1746	count = (/ 1 /) )
1747	CALL handle_netcdf_error( 'netcdf', 126 )
1748
1749	!
1750	!-- Write data for x (shifted by +dx/2) and xu axis
1751	ALLOCATE( netcdf_data(0:nx+1) )
1752
1753	DO i = 0, nx+1
1754	netcdf_data(i) = ( i + 0.5 ) * dx
1755	ENDDO
1756
1757	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_x_xy(av), netcdf_data, &
1758	start = (/ 1 /), count = (/ nx+2 /) )
1759	CALL handle_netcdf_error( 'netcdf', 127 )
1760
1761	DO i = 0, nx+1
1762	netcdf_data(i) = i * dx
1763	ENDDO
1764
1765	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_xu_xy(av), &
1766	netcdf_data, start = (/ 1 /), &
1767	count = (/ nx+2 /) )
1768	CALL handle_netcdf_error( 'netcdf', 367 )
1769
1770	DEALLOCATE( netcdf_data )
1771
1772	!
1773	!-- Write data for y (shifted by +dy/2) and yv axis
1774	ALLOCATE( netcdf_data(0:ny+1) )
1775
1776	DO i = 0, ny+1
1777	netcdf_data(i) = ( i + 0.5 ) * dy
1778	ENDDO
1779
1780	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_y_xy(av), netcdf_data, &
1781	start = (/ 1 /), count = (/ ny+2 /))
1782	CALL handle_netcdf_error( 'netcdf', 128 )
1783
1784	DO i = 0, ny+1
1785	netcdf_data(i) = i * dy
1786	ENDDO
1787
1788	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_yv_xy(av), &
1789	netcdf_data, start = (/ 1 /), &
1790	count = (/ ny+2 /))
1791	CALL handle_netcdf_error( 'netcdf', 368 )
1792
1793	DEALLOCATE( netcdf_data )
1794
1795	!
1796	!-- In case of non-flat topography write height information
1797	IF ( TRIM( topography ) /= 'flat' ) THEN
1798
1799	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av), &
1800	zu_s_inner(0:nx+1,0:ny+1), &
1801	start = (/ 1, 1 /), &
1802	count = (/ nx+2, ny+2 /) )
1803	CALL handle_netcdf_error( 'netcdf', 427 )
1804
1805	nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av), &
1806	zw_w_inner(0:nx+1,0:ny+1), &
1807	start = (/ 1, 1 /), &
1808	count = (/ nx+2, ny+2 /) )
1809	CALL handle_netcdf_error( 'netcdf', 428 )
1810
1811	ENDIF
1812
1813
1814	CASE ( 'xy_ext' )
1815
1816	!
1817	!-- Get the list of variables and compare with the actual run.
1818	!-- First var_list_old has to be reset, since GET_ATT does not assign
1819	!-- trailing blanks.
1820	var_list_old = ' '
1821	nc_stat = NF90_GET_ATT( id_set_xy(av), NF90_GLOBAL, 'VAR_LIST', &
1822	var_list_old )
1823	CALL handle_netcdf_error( 'netcdf', 129 )
1824
1825	var_list = ';'
1826	i = 1
1827	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
1828	IF ( INDEX( do2d(av,i), 'xy' ) /= 0 ) THEN
1829	netcdf_var_name = do2d(av,i)
1830	CALL clean_netcdf_varname( netcdf_var_name )
1831	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
1832	ENDIF
1833	i = i + 1
1834	ENDDO
1835
1836	IF ( av == 0 ) THEN
1837	var = '(xy)'
1838	ELSE
1839	var = '(xy_av)'
1840	ENDIF
1841
1842	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
1843	message_string = 'NetCDF file for cross-sections ' // &
1844	TRIM( var ) // ' from previous run found,' // &
1845	'& but this file cannot be extended due to' // &
1846	' variable mismatch.' // &
1847	'&New file is created instead.'
1848	CALL message( 'define_netcdf_header', 'PA0249', 0, 1, 0, 6, 0 )
1849	extend = .FALSE.
1850	RETURN
1851	ENDIF
1852
1853	!
1854	!-- Calculate the number of current sections
1855	ns = 1
1856	DO WHILE ( section(ns,1) /= -9999 .AND. ns <= 100 )
1857	ns = ns + 1
1858	ENDDO
1859	ns = ns - 1
1860
1861	!
1862	!-- Get and compare the number of horizontal cross sections
1863	nc_stat = NF90_INQ_VARID( id_set_xy(av), 'zu_xy', id_var_zu_xy(av) )
1864	CALL handle_netcdf_error( 'netcdf', 130 )
1865
1866	nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_zu_xy(av), &
1867	dimids = id_dim_zu_xy_old )
1868	CALL handle_netcdf_error( 'netcdf', 131 )
1869	id_dim_zu_xy(av) = id_dim_zu_xy_old(1)
1870
1871	nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_zu_xy(av), &
1872	len = ns_old )
1873	CALL handle_netcdf_error( 'netcdf', 132 )
1874
1875	IF ( ns /= ns_old ) THEN
1876	message_string = 'NetCDF file for cross-sections ' // &
1877	TRIM( var ) // ' from previous run found,' // &
1878	'&but this file cannot be extended due to' // &
1879	' mismatch in number of' // &
1880	'&cross sections.' // &
1881	'&New file is created instead.'
1882	CALL message( 'define_netcdf_header', 'PA0250', 0, 1, 0, 6, 0 )
1883	extend = .FALSE.
1884	RETURN
1885	ENDIF
1886
1887	!
1888	!-- Get and compare the heights of the cross sections
1889	ALLOCATE( netcdf_data(1:ns_old) )
1890
1891	nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_zu_xy(av), netcdf_data )
1892	CALL handle_netcdf_error( 'netcdf', 133 )
1893
1894	DO i = 1, ns
1895	IF ( section(i,1) /= -1 ) THEN
1896	IF ( zu(section(i,1)) /= netcdf_data(i) ) THEN
1897	message_string = 'NetCDF file for cross-sections ' // &
1898	TRIM( var ) // ' from previous run found,' // &
1899	'&but this file cannot be extended' // &
1900	' due to mismatch in cross' // &
1901	'&section levels.' // &
1902	'&New file is created instead.'
1903	CALL message( 'define_netcdf_header', 'PA0251', &
1904	0, 1, 0, 6, 0 )
1905	extend = .FALSE.
1906	RETURN
1907	ENDIF
1908	ELSE
1909	IF ( -1.0 /= netcdf_data(i) ) THEN
1910	message_string = 'NetCDF file for cross-sections ' // &
1911	TRIM( var ) // ' from previous run found,' // &
1912	'&but this file cannot be extended' // &
1913	' due to mismatch in cross' // &
1914	'&section levels.' // &
1915	'&New file is created instead.'
1916	CALL message( 'define_netcdf_header', 'PA0251',&
1917	0, 1, 0, 6, 0 )
1918	extend = .FALSE.
1919	RETURN
1920	ENDIF
1921	ENDIF
1922	ENDDO
1923
1924	DEALLOCATE( netcdf_data )
1925
1926	!
1927	!-- Get the id of the time coordinate (unlimited coordinate) and its
1928	!-- last index on the file. The next time level is do2d..count+1.
1929	!-- The current time must be larger than the last output time
1930	!-- on the file.
1931	nc_stat = NF90_INQ_VARID( id_set_xy(av), 'time', id_var_time_xy(av) )
1932	CALL handle_netcdf_error( 'netcdf', 134 )
1933
1934	nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_time_xy(av), &
1935	dimids = id_dim_time_old )
1936	CALL handle_netcdf_error( 'netcdf', 135 )
1937	id_dim_time_xy(av) = id_dim_time_old(1)
1938
1939	nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_time_xy(av), &
1940	len = do2d_xy_time_count(av) )
1941	CALL handle_netcdf_error( 'netcdf', 136 )
1942
1943	nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_time_xy(av), &
1944	last_time_coordinate, &
1945	start = (/ do2d_xy_time_count(av) /), &
1946	count = (/ 1 /) )
1947	CALL handle_netcdf_error( 'netcdf', 137 )
1948
1949	IF ( last_time_coordinate(1) >= simulated_time ) THEN
1950	message_string = 'NetCDF file for cross sections ' // &
1951	TRIM( var ) // ' from previous run found,' // &
1952	'&but this file cannot be extended becaus' // &
1953	'e the current output time' // &
1954	'&is less or equal than the last output t' // &
1955	'ime on this file.' // &
1956	'&New file is created instead.'
1957	CALL message( 'define_netcdf_header', 'PA0252', 0, 1, 0, 6, 0 )
1958	do2d_xy_time_count(av) = 0
1959	extend = .FALSE.
1960	RETURN
1961	ENDIF
1962
1963	!
1964	!-- Dataset seems to be extendable.
1965	!-- Now get the variable ids.
1966	i = 1
1967	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
1968	IF ( INDEX( do2d(av,i), 'xy' ) /= 0 ) THEN
1969	netcdf_var_name = do2d(av,i)
1970	CALL clean_netcdf_varname( netcdf_var_name )
1971	nc_stat = NF90_INQ_VARID( id_set_xy(av), netcdf_var_name, &
1972	id_var_do2d(av,i) )
1973	CALL handle_netcdf_error( 'netcdf', 138 )
1974	#if defined( __netcdf4 )
1975	!
1976	!-- Set collective io operations for parallel io
1977	IF ( netcdf_data_format > 2 ) THEN
1978	nc_stat = NF90_VAR_PAR_ACCESS( id_set_xy(av), &
1979	id_var_do2d(av,i), &
1980	NF90_COLLECTIVE )
1981	CALL handle_netcdf_error( 'netcdf', 454 )
1982	ENDIF
1983	#endif
1984	ENDIF
1985	i = i + 1
1986	ENDDO
1987
1988	!
1989	!-- Update the title attribute on file
1990	!-- In order to avoid 'data mode' errors if updated attributes are larger
1991	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
1992	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
1993	!-- performance loss due to data copying; an alternative strategy would be
1994	!-- to ensure equal attribute size in a job chain. Maybe revise later.
1995	IF ( av == 0 ) THEN
1996	time_average_text = ' '
1997	ELSE
1998	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
1999	averaging_interval
2000	ENDIF
2001	nc_stat = NF90_REDEF( id_set_xy(av) )
2002	CALL handle_netcdf_error( 'netcdf', 431 )
2003	nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'title', &
2004	TRIM( run_description_header ) // &
2005	TRIM( time_average_text ) )
2006	CALL handle_netcdf_error( 'netcdf', 139 )
2007	nc_stat = NF90_ENDDEF( id_set_xy(av) )
2008	CALL handle_netcdf_error( 'netcdf', 432 )
2009	message_string = 'NetCDF file for cross-sections ' // &
2010	TRIM( var ) // ' from previous run found.' // &
2011	'&This file will be extended.'
2012	CALL message( 'define_netcdf_header', 'PA0253', 0, 0, 0, 6, 0 )
2013
2014
2015	CASE ( 'xz_new' )
2016
2017	!
2018	!-- Define some global attributes of the dataset
2019	nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'Conventions', &
2020	'COARDS' )
2021	CALL handle_netcdf_error( 'netcdf', 140 )
2022
2023	IF ( av == 0 ) THEN
2024	time_average_text = ' '
2025	ELSE
2026	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
2027	averaging_interval
2028	ENDIF
2029	nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'title', &
2030	TRIM( run_description_header ) // &
2031	TRIM( time_average_text ) )
2032	CALL handle_netcdf_error( 'netcdf', 141 )
2033	IF ( av == 1 ) THEN
2034	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval
2035	nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'time_avg', &
2036	TRIM( time_average_text ) )
2037	CALL handle_netcdf_error( 'netcdf', 141 )
2038	ENDIF
2039
2040	!
2041	!-- Define time coordinate for xz sections (unlimited dimension)
2042	nc_stat = NF90_DEF_DIM( id_set_xz(av), 'time', NF90_UNLIMITED, &
2043	id_dim_time_xz(av) )
2044	CALL handle_netcdf_error( 'netcdf', 142 )
2045
2046	nc_stat = NF90_DEF_VAR( id_set_xz(av), 'time', NF90_DOUBLE, &
2047	id_dim_time_xz(av), id_var_time_xz(av) )
2048	CALL handle_netcdf_error( 'netcdf', 143 )
2049
2050	nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_time_xz(av), 'units', &
2051	'seconds')
2052	CALL handle_netcdf_error( 'netcdf', 144 )
2053
2054	!
2055	!-- Define the spatial dimensions and coordinates for xz-sections.
2056	!-- First, determine the number of vertical sections to be written.
2057	IF ( section(1,2) == -9999 ) THEN
2058	RETURN
2059	ELSE
2060	ns = 1
2061	DO WHILE ( section(ns,2) /= -9999 .AND. ns <= 100 )
2062	ns = ns + 1
2063	ENDDO
2064	ns = ns - 1
2065	ENDIF
2066
2067	!
2068	!-- Define y-axis (for scalar position)
2069	nc_stat = NF90_DEF_DIM( id_set_xz(av), 'y_xz', ns, id_dim_y_xz(av) )
2070	CALL handle_netcdf_error( 'netcdf', 145 )
2071
2072	nc_stat = NF90_DEF_VAR( id_set_xz(av), 'y_xz', NF90_DOUBLE, &
2073	id_dim_y_xz(av), id_var_y_xz(av) )
2074	CALL handle_netcdf_error( 'netcdf', 146 )
2075
2076	nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_y_xz(av), 'units', &
2077	'meters' )
2078	CALL handle_netcdf_error( 'netcdf', 147 )
2079
2080	!
2081	!-- Define y-axis (for v position)
2082	nc_stat = NF90_DEF_DIM( id_set_xz(av), 'yv_xz', ns, id_dim_yv_xz(av) )
2083	CALL handle_netcdf_error( 'netcdf', 369 )
2084
2085	nc_stat = NF90_DEF_VAR( id_set_xz(av), 'yv_xz', NF90_DOUBLE, &
2086	id_dim_yv_xz(av), id_var_yv_xz(av) )
2087	CALL handle_netcdf_error( 'netcdf', 370 )
2088
2089	nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_yv_xz(av), 'units', &
2090	'meters' )
2091	CALL handle_netcdf_error( 'netcdf', 371 )
2092
2093	!
2094	!-- Define a variable to store the layer indices of the vertical cross
2095	!-- sections
2096	nc_stat = NF90_DEF_VAR( id_set_xz(av), 'ind_y_xz', NF90_DOUBLE, &
2097	id_dim_y_xz(av), id_var_ind_y_xz(av) )
2098	CALL handle_netcdf_error( 'netcdf', 148 )
2099
2100	nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_ind_y_xz(av), 'units', &
2101	'gridpoints')
2102	CALL handle_netcdf_error( 'netcdf', 149 )
2103
2104	!
2105	!-- Define x-axis (for scalar position)
2106	nc_stat = NF90_DEF_DIM( id_set_xz(av), 'x', nx+2, id_dim_x_xz(av) )
2107	CALL handle_netcdf_error( 'netcdf', 150 )
2108
2109	nc_stat = NF90_DEF_VAR( id_set_xz(av), 'x', NF90_DOUBLE, &
2110	id_dim_x_xz(av), id_var_x_xz(av) )
2111	CALL handle_netcdf_error( 'netcdf', 151 )
2112
2113	nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_x_xz(av), 'units', &
2114	'meters' )
2115	CALL handle_netcdf_error( 'netcdf', 152 )
2116
2117	!
2118	!-- Define x-axis (for u position)
2119	nc_stat = NF90_DEF_DIM( id_set_xz(av), 'xu', nx+2, id_dim_xu_xz(av) )
2120	CALL handle_netcdf_error( 'netcdf', 372 )
2121
2122	nc_stat = NF90_DEF_VAR( id_set_xz(av), 'xu', NF90_DOUBLE, &
2123	id_dim_xu_xz(av), id_var_xu_xz(av) )
2124	CALL handle_netcdf_error( 'netcdf', 373 )
2125
2126	nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_xu_xz(av), 'units', &
2127	'meters' )
2128	CALL handle_netcdf_error( 'netcdf', 374 )
2129
2130	!
2131	!-- Define the two z-axes (zu and zw)
2132	nc_stat = NF90_DEF_DIM( id_set_xz(av), 'zu', nz+2, id_dim_zu_xz(av) )
2133	CALL handle_netcdf_error( 'netcdf', 153 )
2134
2135	nc_stat = NF90_DEF_VAR( id_set_xz(av), 'zu', NF90_DOUBLE, &
2136	id_dim_zu_xz(av), id_var_zu_xz(av) )
2137	CALL handle_netcdf_error( 'netcdf', 154 )
2138
2139	nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_zu_xz(av), 'units', &
2140	'meters' )
2141	CALL handle_netcdf_error( 'netcdf', 155 )
2142
2143	nc_stat = NF90_DEF_DIM( id_set_xz(av), 'zw', nz+2, id_dim_zw_xz(av) )
2144	CALL handle_netcdf_error( 'netcdf', 156 )
2145
2146	nc_stat = NF90_DEF_VAR( id_set_xz(av), 'zw', NF90_DOUBLE, &
2147	id_dim_zw_xz(av), id_var_zw_xz(av) )
2148	CALL handle_netcdf_error( 'netcdf', 157 )
2149
2150	nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_zw_xz(av), 'units', &
2151	'meters' )
2152	CALL handle_netcdf_error( 'netcdf', 158 )
2153
2154
2155	!
2156	!-- Define the variables
2157	var_list = ';'
2158	i = 1
2159
2160	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2161
2162	IF ( INDEX( do2d(av,i), 'xz' ) /= 0 ) THEN
2163
2164	!
2165	!-- Check for the grid
2166	found = .TRUE.
2167	SELECT CASE ( do2d(av,i) )
2168	!
2169	!-- Most variables are defined on the zu grid
2170	CASE ( 'e_xz', 'p_xz', 'pc_xz', 'pr_xz', 'pt_xz', 'q_xz', &
2171	'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qv_xz', &
2172	'rho_xz', 's_xz', 'sa_xz', 'vpt_xz' )
2173
2174	grid_x = 'x'
2175	grid_y = 'y'
2176	grid_z = 'zu'
2177	!
2178	!-- u grid
2179	CASE ( 'u_xz' )
2180
2181	grid_x = 'xu'
2182	grid_y = 'y'
2183	grid_z = 'zu'
2184	!
2185	!-- v grid
2186	CASE ( 'v_xz' )
2187
2188	grid_x = 'x'
2189	grid_y = 'yv'
2190	grid_z = 'zu'
2191	!
2192	!-- w grid
2193	CASE ( 'w_xz' )
2194
2195	grid_x = 'x'
2196	grid_y = 'y'
2197	grid_z = 'zw'
2198
2199	CASE DEFAULT
2200	!
2201	!-- Check for user-defined quantities
2202	CALL user_define_netcdf_grid( do2d(av,i), found, &
2203	grid_x, grid_y, grid_z )
2204
2205	END SELECT
2206
2207	!
2208	!-- Select the respective dimension ids
2209	IF ( grid_x == 'x' ) THEN
2210	id_x = id_dim_x_xz(av)
2211	ELSEIF ( grid_x == 'xu' ) THEN
2212	id_x = id_dim_xu_xz(av)
2213	ENDIF
2214
2215	IF ( grid_y == 'y' ) THEN
2216	id_y = id_dim_y_xz(av)
2217	ELSEIF ( grid_y == 'yv' ) THEN
2218	id_y = id_dim_yv_xz(av)
2219	ENDIF
2220
2221	IF ( grid_z == 'zu' ) THEN
2222	id_z = id_dim_zu_xz(av)
2223	ELSEIF ( grid_z == 'zw' ) THEN
2224	id_z = id_dim_zw_xz(av)
2225	ENDIF
2226
2227	!
2228	!-- Define the grid
2229	nc_stat = NF90_DEF_VAR( id_set_xz(av), do2d(av,i), &
2230	nc_precision(2), &
2231	(/ id_x, id_y, id_z, id_dim_time_xz(av) /), &
2232	id_var_do2d(av,i) )
2233
2234	IF ( .NOT. found ) THEN
2235	WRITE ( message_string, * ) 'no grid defined for', &
2236	' variable ', TRIM( do2d(av,i) )
2237	CALL message( 'define_netcdf_header', 'PA0244',&
2238	0, 1, 0, 6, 0 )
2239	ENDIF
2240
2241	var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
2242
2243	CALL handle_netcdf_error( 'netcdf', 159 )
2244	!
2245	!-- Store the 'real' name of the variable (with *, for example)
2246	!-- in the long_name attribute. This is evaluated by Ferret,
2247	!-- for example.
2248	nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_do2d(av,i), &
2249	'long_name', do2d(av,i) )
2250	CALL handle_netcdf_error( 'netcdf', 160 )
2251	!
2252	!-- Define the variable's unit
2253	nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_do2d(av,i), &
2254	'units', TRIM( do2d_unit(av,i) ) )
2255	CALL handle_netcdf_error( 'netcdf', 355 )
2256	#if defined( __netcdf4 )
2257	!
2258	!-- Set independent io operations for parallel io. Collective io
2259	!-- is only allowed in case of a 1d-decomposition along x, because
2260	!-- otherwise, not all PEs have output data.
2261	IF ( netcdf_data_format > 2 ) THEN
2262	IF ( npey == 1 ) THEN
2263	nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av), &
2264	id_var_do2d(av,i), &
2265	NF90_COLLECTIVE )
2266	ELSE
2267	!
2268	!-- ATTENTION: Due to a probable bug in the NetCDF4
2269	!-- installation, independet output does not work
2270	!-- A workaround is used in data_output_2d on those
2271	!-- PEs having no data
2272	nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av), &
2273	id_var_do2d(av,i), &
2274	NF90_COLLECTIVE )
2275	! nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av), &
2276	! id_var_do2d(av,i), &
2277	! NF90_INDEPENDENT )
2278	ENDIF
2279	CALL handle_netcdf_error( 'netcdf', 449 )
2280	ENDIF
2281	#endif
2282	ENDIF
2283
2284	i = i + 1
2285
2286	ENDDO
2287
2288	!
2289	!-- No arrays to output. Close the netcdf file and return.
2290	IF ( i == 1 ) RETURN
2291
2292	!
2293	!-- Write the list of variables as global attribute (this is used by
2294	!-- restart runs and by combine_plot_fields)
2295	nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'VAR_LIST', &
2296	var_list )
2297	CALL handle_netcdf_error( 'netcdf', 161 )
2298
2299	!
2300	!-- Leave NetCDF define mode
2301	nc_stat = NF90_ENDDEF( id_set_xz(av) )
2302	CALL handle_netcdf_error( 'netcdf', 162 )
2303
2304	!
2305	!-- Write axis data: y_xz, x, zu, zw
2306	ALLOCATE( netcdf_data(1:ns) )
2307
2308	!
2309	!-- Write y_xz data (shifted by +dy/2)
2310	DO i = 1, ns
2311	IF( section(i,2) == -1 ) THEN
2312	netcdf_data(i) = -1.0 ! section averaged along y
2313	ELSE
2314	netcdf_data(i) = ( section(i,2) + 0.5 ) * dy
2315	ENDIF
2316	ENDDO
2317	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_y_xz(av), netcdf_data, &
2318	start = (/ 1 /), count = (/ ns /) )
2319	CALL handle_netcdf_error( 'netcdf', 163 )
2320
2321	!
2322	!-- Write yv_xz data
2323	DO i = 1, ns
2324	IF( section(i,2) == -1 ) THEN
2325	netcdf_data(i) = -1.0 ! section averaged along y
2326	ELSE
2327	netcdf_data(i) = section(i,2) * dy
2328	ENDIF
2329	ENDDO
2330	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_yv_xz(av), &
2331	netcdf_data, start = (/ 1 /), &
2332	count = (/ ns /) )
2333	CALL handle_netcdf_error( 'netcdf', 375 )
2334
2335	!
2336	!-- Write gridpoint number data
2337	netcdf_data(1:ns) = section(1:ns,2)
2338	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_ind_y_xz(av), &
2339	netcdf_data, start = (/ 1 /), &
2340	count = (/ ns /) )
2341	CALL handle_netcdf_error( 'netcdf', 164 )
2342
2343
2344	DEALLOCATE( netcdf_data )
2345
2346	!
2347	!-- Write data for x (shifted by +dx/2) and xu axis
2348	ALLOCATE( netcdf_data(0:nx+1) )
2349
2350	DO i = 0, nx+1
2351	netcdf_data(i) = ( i + 0.5 ) * dx
2352	ENDDO
2353
2354	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_x_xz(av), netcdf_data, &
2355	start = (/ 1 /), count = (/ nx+2 /) )
2356	CALL handle_netcdf_error( 'netcdf', 165 )
2357
2358	DO i = 0, nx+1
2359	netcdf_data(i) = i * dx
2360	ENDDO
2361
2362	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_xu_xz(av), &
2363	netcdf_data, start = (/ 1 /), &
2364	count = (/ nx+2 /) )
2365	CALL handle_netcdf_error( 'netcdf', 377 )
2366
2367	DEALLOCATE( netcdf_data )
2368
2369	!
2370	!-- Write zu and zw data (vertical axes)
2371	ALLOCATE( netcdf_data(0:nz+1) )
2372
2373	netcdf_data(0:nz+1) = zu(nzb:nzt+1)
2374	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zu_xz(av), &
2375	netcdf_data, start = (/ 1 /), &
2376	count = (/ nz+2 /) )
2377	CALL handle_netcdf_error( 'netcdf', 166 )
2378
2379	netcdf_data(0:nz+1) = zw(nzb:nzt+1)
2380	nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zw_xz(av), &
2381	netcdf_data, start = (/ 1 /), &
2382	count = (/ nz+2 /) )
2383	CALL handle_netcdf_error( 'netcdf', 167 )
2384
2385	DEALLOCATE( netcdf_data )
2386
2387
2388	CASE ( 'xz_ext' )
2389
2390	!
2391	!-- Get the list of variables and compare with the actual run.
2392	!-- First var_list_old has to be reset, since GET_ATT does not assign
2393	!-- trailing blanks.
2394	var_list_old = ' '
2395	nc_stat = NF90_GET_ATT( id_set_xz(av), NF90_GLOBAL, 'VAR_LIST', &
2396	var_list_old )
2397	CALL handle_netcdf_error( 'netcdf', 168 )
2398
2399	var_list = ';'
2400	i = 1
2401	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2402	IF ( INDEX( do2d(av,i), 'xz' ) /= 0 ) THEN
2403	netcdf_var_name = do2d(av,i)
2404	CALL clean_netcdf_varname( netcdf_var_name )
2405	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
2406	ENDIF
2407	i = i + 1
2408	ENDDO
2409
2410	IF ( av == 0 ) THEN
2411	var = '(xz)'
2412	ELSE
2413	var = '(xz_av)'
2414	ENDIF
2415
2416	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
2417	message_string = 'NetCDF file for cross-sections ' // &
2418	TRIM( var ) // ' from previous run found,' // &
2419	'& but this file cannot be extended due to' // &
2420	' variable mismatch.' // &
2421	'&New file is created instead.'
2422	CALL message( 'define_netcdf_header', 'PA0249', 0, 1, 0, 6, 0 )
2423	extend = .FALSE.
2424	RETURN
2425	ENDIF
2426
2427	!
2428	!-- Calculate the number of current sections
2429	ns = 1
2430	DO WHILE ( section(ns,2) /= -9999 .AND. ns <= 100 )
2431	ns = ns + 1
2432	ENDDO
2433	ns = ns - 1
2434
2435	!
2436	!-- Get and compare the number of vertical cross sections
2437	nc_stat = NF90_INQ_VARID( id_set_xz(av), 'y_xz', id_var_y_xz(av) )
2438	CALL handle_netcdf_error( 'netcdf', 169 )
2439
2440	nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_y_xz(av), &
2441	dimids = id_dim_y_xz_old )
2442	CALL handle_netcdf_error( 'netcdf', 170 )
2443	id_dim_y_xz(av) = id_dim_y_xz_old(1)
2444
2445	nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_y_xz(av), &
2446	len = ns_old )
2447	CALL handle_netcdf_error( 'netcdf', 171 )
2448
2449	IF ( ns /= ns_old ) THEN
2450	message_string = 'NetCDF file for cross-sections ' // &
2451	TRIM( var ) // ' from previous run found,' // &
2452	'&but this file cannot be extended due to' // &
2453	' mismatch in number of' // &
2454	'&cross sections.' // &
2455	'&New file is created instead.'
2456	CALL message( 'define_netcdf_header', 'PA0250', 0, 1, 0, 6, 0 )
2457	extend = .FALSE.
2458	RETURN
2459	ENDIF
2460
2461	!
2462	!-- Get and compare the heights of the cross sections
2463	ALLOCATE( netcdf_data(1:ns_old) )
2464
2465	nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_y_xz(av), netcdf_data )
2466	CALL handle_netcdf_error( 'netcdf', 172 )
2467
2468	DO i = 1, ns
2469	IF ( section(i,2) /= -1 ) THEN
2470	IF ( ( section(i,2) * dy ) /= netcdf_data(i) ) THEN
2471	message_string = 'NetCDF file for cross-sections ' // &
2472	TRIM( var ) // ' from previous run found,' // &
2473	'&but this file cannot be extended' // &
2474	' due to mismatch in cross' // &
2475	'&section levels.' // &
2476	'&New file is created instead.'
2477	CALL message( 'define_netcdf_header', 'PA0251', &
2478	0, 1, 0, 6, 0 )
2479	extend = .FALSE.
2480	RETURN
2481	ENDIF
2482	ELSE
2483	IF ( -1.0 /= netcdf_data(i) ) THEN
2484	message_string = 'NetCDF file for cross-sections ' // &
2485	TRIM( var ) // ' from previous run found,' // &
2486	'&but this file cannot be extended' // &
2487	' due to mismatch in cross' // &
2488	'&section levels.' // &
2489	'&New file is created instead.'
2490	CALL message( 'define_netcdf_header', 'PA0251', &
2491	0, 1, 0, 6, 0 )
2492	extend = .FALSE.
2493	RETURN
2494	ENDIF
2495	ENDIF
2496	ENDDO
2497
2498	DEALLOCATE( netcdf_data )
2499
2500	!
2501	!-- Get the id of the time coordinate (unlimited coordinate) and its
2502	!-- last index on the file. The next time level is do2d..count+1.
2503	!-- The current time must be larger than the last output time
2504	!-- on the file.
2505	nc_stat = NF90_INQ_VARID( id_set_xz(av), 'time', id_var_time_xz(av) )
2506	CALL handle_netcdf_error( 'netcdf', 173 )
2507
2508	nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_time_xz(av), &
2509	dimids = id_dim_time_old )
2510	CALL handle_netcdf_error( 'netcdf', 174 )
2511	id_dim_time_xz(av) = id_dim_time_old(1)
2512
2513	nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_time_xz(av), &
2514	len = do2d_xz_time_count(av) )
2515	CALL handle_netcdf_error( 'netcdf', 175 )
2516
2517	nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_time_xz(av), &
2518	last_time_coordinate, &
2519	start = (/ do2d_xz_time_count(av) /), &
2520	count = (/ 1 /) )
2521	CALL handle_netcdf_error( 'netcdf', 176 )
2522
2523	IF ( last_time_coordinate(1) >= simulated_time ) THEN
2524	message_string = 'NetCDF file for cross sections ' // &
2525	TRIM( var ) // ' from previous run found,' // &
2526	'&but this file cannot be extended becaus' // &
2527	'e the current output time' // &
2528	'&is less or equal than the last output t' // &
2529	'ime on this file.' // &
2530	'&New file is created instead.'
2531	CALL message( 'define_netcdf_header', 'PA0252', 0, 1, 0, 6, 0 )
2532	do2d_xz_time_count(av) = 0
2533	extend = .FALSE.
2534	RETURN
2535	ENDIF
2536
2537	!
2538	!-- Dataset seems to be extendable.
2539	!-- Now get the variable ids.
2540	i = 1
2541	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2542	IF ( INDEX( do2d(av,i), 'xz' ) /= 0 ) THEN
2543	netcdf_var_name = do2d(av,i)
2544	CALL clean_netcdf_varname( netcdf_var_name )
2545	nc_stat = NF90_INQ_VARID( id_set_xz(av), netcdf_var_name, &
2546	id_var_do2d(av,i) )
2547	CALL handle_netcdf_error( 'netcdf', 177 )
2548	#if defined( __netcdf4 )
2549	!
2550	!-- Set independent io operations for parallel io. Collective io
2551	!-- is only allowed in case of a 1d-decomposition along x, because
2552	!-- otherwise, not all PEs have output data.
2553	IF ( netcdf_data_format > 2 ) THEN
2554	IF ( npey == 1 ) THEN
2555	nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av), &
2556	id_var_do2d(av,i), &
2557	NF90_COLLECTIVE )
2558	ELSE
2559	!
2560	!-- ATTENTION: Due to a probable bug in the NetCDF4
2561	!-- installation, independet output does not work
2562	!-- A workaround is used in data_output_2d on those
2563	!-- PEs having no data
2564	nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av), &
2565	id_var_do2d(av,i), &
2566	NF90_COLLECTIVE )
2567	! nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av), &
2568	! id_var_do2d(av,i), &
2569	! NF90_INDEPENDENT )
2570	ENDIF
2571	CALL handle_netcdf_error( 'netcdf', 455 )
2572	ENDIF
2573	#endif
2574	ENDIF
2575	i = i + 1
2576	ENDDO
2577
2578	!
2579	!-- Update the title attribute on file
2580	!-- In order to avoid 'data mode' errors if updated attributes are larger
2581	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
2582	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
2583	!-- performance loss due to data copying; an alternative strategy would be
2584	!-- to ensure equal attribute size in a job chain. Maybe revise later.
2585	IF ( av == 0 ) THEN
2586	time_average_text = ' '
2587	ELSE
2588	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
2589	averaging_interval
2590	ENDIF
2591	nc_stat = NF90_REDEF( id_set_xz(av) )
2592	CALL handle_netcdf_error( 'netcdf', 433 )
2593	nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'title', &
2594	TRIM( run_description_header ) // &
2595	TRIM( time_average_text ) )
2596	CALL handle_netcdf_error( 'netcdf', 178 )
2597	nc_stat = NF90_ENDDEF( id_set_xz(av) )
2598	CALL handle_netcdf_error( 'netcdf', 434 )
2599	message_string = 'NetCDF file for cross-sections ' // &
2600	TRIM( var ) // ' from previous run found.' // &
2601	'&This file will be extended.'
2602	CALL message( 'define_netcdf_header', 'PA0253', 0, 0, 0, 6, 0 )
2603
2604
2605	CASE ( 'yz_new' )
2606
2607	!
2608	!-- Define some global attributes of the dataset
2609	nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'Conventions', &
2610	'COARDS' )
2611	CALL handle_netcdf_error( 'netcdf', 179 )
2612
2613	IF ( av == 0 ) THEN
2614	time_average_text = ' '
2615	ELSE
2616	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
2617	averaging_interval
2618	ENDIF
2619	nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'title', &
2620	TRIM( run_description_header ) // &
2621	TRIM( time_average_text ) )
2622	CALL handle_netcdf_error( 'netcdf', 180 )
2623	IF ( av == 1 ) THEN
2624	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval
2625	nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'time_avg', &
2626	TRIM( time_average_text ) )
2627	CALL handle_netcdf_error( 'netcdf', 180 )
2628	ENDIF
2629
2630	!
2631	!-- Define time coordinate for yz sections (unlimited dimension)
2632	nc_stat = NF90_DEF_DIM( id_set_yz(av), 'time', NF90_UNLIMITED, &
2633	id_dim_time_yz(av) )
2634	CALL handle_netcdf_error( 'netcdf', 181 )
2635
2636	nc_stat = NF90_DEF_VAR( id_set_yz(av), 'time', NF90_DOUBLE, &
2637	id_dim_time_yz(av), id_var_time_yz(av) )
2638	CALL handle_netcdf_error( 'netcdf', 182 )
2639
2640	nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_time_yz(av), 'units', &
2641	'seconds')
2642	CALL handle_netcdf_error( 'netcdf', 183 )
2643
2644	!
2645	!-- Define the spatial dimensions and coordinates for yz-sections.
2646	!-- First, determine the number of vertical sections to be written.
2647	IF ( section(1,3) == -9999 ) THEN
2648	RETURN
2649	ELSE
2650	ns = 1
2651	DO WHILE ( section(ns,3) /= -9999 .AND. ns <= 100 )
2652	ns = ns + 1
2653	ENDDO
2654	ns = ns - 1
2655	ENDIF
2656
2657	!
2658	!-- Define x axis (for scalar position)
2659	nc_stat = NF90_DEF_DIM( id_set_yz(av), 'x_yz', ns, id_dim_x_yz(av) )
2660	CALL handle_netcdf_error( 'netcdf', 184 )
2661
2662	nc_stat = NF90_DEF_VAR( id_set_yz(av), 'x_yz', NF90_DOUBLE, &
2663	id_dim_x_yz(av), id_var_x_yz(av) )
2664	CALL handle_netcdf_error( 'netcdf', 185 )
2665
2666	nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_x_yz(av), 'units', &
2667	'meters' )
2668	CALL handle_netcdf_error( 'netcdf', 186 )
2669
2670	!
2671	!-- Define x axis (for u position)
2672	nc_stat = NF90_DEF_DIM( id_set_yz(av), 'xu_yz', ns, id_dim_xu_yz(av) )
2673	CALL handle_netcdf_error( 'netcdf', 377 )
2674
2675	nc_stat = NF90_DEF_VAR( id_set_yz(av), 'xu_yz', NF90_DOUBLE, &
2676	id_dim_xu_yz(av), id_var_xu_yz(av) )
2677	CALL handle_netcdf_error( 'netcdf', 378 )
2678
2679	nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_xu_yz(av), 'units', &
2680	'meters' )
2681	CALL handle_netcdf_error( 'netcdf', 379 )
2682
2683	!
2684	!-- Define a variable to store the layer indices of the vertical cross
2685	!-- sections
2686	nc_stat = NF90_DEF_VAR( id_set_yz(av), 'ind_x_yz', NF90_DOUBLE, &
2687	id_dim_x_yz(av), id_var_ind_x_yz(av) )
2688	CALL handle_netcdf_error( 'netcdf', 187 )
2689
2690	nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_ind_x_yz(av), 'units', &
2691	'gridpoints')
2692	CALL handle_netcdf_error( 'netcdf', 188 )
2693
2694	!
2695	!-- Define y-axis (for scalar position)
2696	nc_stat = NF90_DEF_DIM( id_set_yz(av), 'y', ny+2, id_dim_y_yz(av) )
2697	CALL handle_netcdf_error( 'netcdf', 189 )
2698
2699	nc_stat = NF90_DEF_VAR( id_set_yz(av), 'y', NF90_DOUBLE, &
2700	id_dim_y_yz(av), id_var_y_yz(av) )
2701	CALL handle_netcdf_error( 'netcdf', 190 )
2702
2703	nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_y_yz(av), 'units', &
2704	'meters' )
2705	CALL handle_netcdf_error( 'netcdf', 191 )
2706
2707	!
2708	!-- Define y-axis (for v position)
2709	nc_stat = NF90_DEF_DIM( id_set_yz(av), 'yv', ny+2, id_dim_yv_yz(av) )
2710	CALL handle_netcdf_error( 'netcdf', 380 )
2711
2712	nc_stat = NF90_DEF_VAR( id_set_yz(av), 'yv', NF90_DOUBLE, &
2713	id_dim_yv_yz(av), id_var_yv_yz(av) )
2714	CALL handle_netcdf_error( 'netcdf', 381 )
2715
2716	nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_yv_yz(av), 'units', &
2717	'meters' )
2718	CALL handle_netcdf_error( 'netcdf', 382 )
2719
2720	!
2721	!-- Define the two z-axes (zu and zw)
2722	nc_stat = NF90_DEF_DIM( id_set_yz(av), 'zu', nz+2, id_dim_zu_yz(av) )
2723	CALL handle_netcdf_error( 'netcdf', 192 )
2724
2725	nc_stat = NF90_DEF_VAR( id_set_yz(av), 'zu', NF90_DOUBLE, &
2726	id_dim_zu_yz(av), id_var_zu_yz(av) )
2727	CALL handle_netcdf_error( 'netcdf', 193 )
2728
2729	nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_zu_yz(av), 'units', &
2730	'meters' )
2731	CALL handle_netcdf_error( 'netcdf', 194 )
2732
2733	nc_stat = NF90_DEF_DIM( id_set_yz(av), 'zw', nz+2, id_dim_zw_yz(av) )
2734	CALL handle_netcdf_error( 'netcdf', 195 )
2735
2736	nc_stat = NF90_DEF_VAR( id_set_yz(av), 'zw', NF90_DOUBLE, &
2737	id_dim_zw_yz(av), id_var_zw_yz(av) )
2738	CALL handle_netcdf_error( 'netcdf', 196 )
2739
2740	nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_zw_yz(av), 'units', &
2741	'meters' )
2742	CALL handle_netcdf_error( 'netcdf', 197 )
2743
2744
2745	!
2746	!-- Define the variables
2747	var_list = ';'
2748	i = 1
2749
2750	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2751
2752	IF ( INDEX( do2d(av,i), 'yz' ) /= 0 ) THEN
2753
2754	!
2755	!-- Check for the grid
2756	found = .TRUE.
2757	SELECT CASE ( do2d(av,i) )
2758	!
2759	!-- Most variables are defined on the zu grid
2760	CASE ( 'e_yz', 'p_yz', 'pc_yz', 'pr_yz', 'pt_yz', 'q_yz', &
2761	'ql_yz', 'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qv_yz', &
2762	'rho_yz', 's_yz', 'sa_yz', 'vpt_yz' )
2763
2764	grid_x = 'x'
2765	grid_y = 'y'
2766	grid_z = 'zu'
2767	!
2768	!-- u grid
2769	CASE ( 'u_yz' )
2770
2771	grid_x = 'xu'
2772	grid_y = 'y'
2773	grid_z = 'zu'
2774	!
2775	!-- v grid
2776	CASE ( 'v_yz' )
2777
2778	grid_x = 'x'
2779	grid_y = 'yv'
2780	grid_z = 'zu'
2781	!
2782	!-- w grid
2783	CASE ( 'w_yz' )
2784
2785	grid_x = 'x'
2786	grid_y = 'y'
2787	grid_z = 'zw'
2788
2789	CASE DEFAULT
2790	!
2791	!-- Check for user-defined quantities
2792	CALL user_define_netcdf_grid( do2d(av,i), found, &
2793	grid_x, grid_y, grid_z )
2794
2795	END SELECT
2796
2797	!
2798	!-- Select the respective dimension ids
2799	IF ( grid_x == 'x' ) THEN
2800	id_x = id_dim_x_yz(av)
2801	ELSEIF ( grid_x == 'xu' ) THEN
2802	id_x = id_dim_xu_yz(av)
2803	ENDIF
2804
2805	IF ( grid_y == 'y' ) THEN
2806	id_y = id_dim_y_yz(av)
2807	ELSEIF ( grid_y == 'yv' ) THEN
2808	id_y = id_dim_yv_yz(av)
2809	ENDIF
2810
2811	IF ( grid_z == 'zu' ) THEN
2812	id_z = id_dim_zu_yz(av)
2813	ELSEIF ( grid_z == 'zw' ) THEN
2814	id_z = id_dim_zw_yz(av)
2815	ENDIF
2816
2817	!
2818	!-- Define the grid
2819	nc_stat = NF90_DEF_VAR( id_set_yz(av), do2d(av,i), &
2820	nc_precision(3), &
2821	(/ id_x, id_y, id_z, id_dim_time_yz(av) /), &
2822	id_var_do2d(av,i) )
2823
2824	IF ( .NOT. found ) THEN
2825	WRITE ( message_string, * ) 'no grid defined for', &
2826	' variable ', TRIM( do2d(av,i) )
2827	CALL message( 'define_netcdf_header', 'PA0244',&
2828	0, 1, 0, 6, 0 )
2829	ENDIF
2830
2831	var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
2832
2833	CALL handle_netcdf_error( 'netcdf', 198 )
2834	!
2835	!-- Store the 'real' name of the variable (with *, for example)
2836	!-- in the long_name attribute. This is evaluated by Ferret,
2837	!-- for example.
2838	nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_do2d(av,i), &
2839	'long_name', do2d(av,i) )
2840	CALL handle_netcdf_error( 'netcdf', 199 )
2841	!
2842	!-- Define the variable's unit
2843	nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_do2d(av,i), &
2844	'units', TRIM( do2d_unit(av,i) ) )
2845	CALL handle_netcdf_error( 'netcdf', 356 )
2846	#if defined( __netcdf4 )
2847	!
2848	!-- Set independent io operations for parallel io. Collective io
2849	!-- is only allowed in case of a 1d-decomposition along y, because
2850	!-- otherwise, not all PEs have output data.
2851	IF ( netcdf_data_format > 2 ) THEN
2852	IF ( npex == 1 ) THEN
2853	nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av), &
2854	id_var_do2d(av,i), &
2855	NF90_COLLECTIVE )
2856	ELSE
2857	!
2858	!-- ATTENTION: Due to a probable bug in the NetCDF4
2859	!-- installation, independet output does not work
2860	!-- A workaround is used in data_output_2d on those
2861	!-- PEs having no data
2862	nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av), &
2863	id_var_do2d(av,i), &
2864	NF90_COLLECTIVE )
2865	! nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av), &
2866	! id_var_do2d(av,i), &
2867	! NF90_INDEPENDENT )
2868	ENDIF
2869	CALL handle_netcdf_error( 'netcdf', 450 )
2870	ENDIF
2871	#endif
2872	ENDIF
2873
2874	i = i + 1
2875
2876	ENDDO
2877
2878	!
2879	!-- No arrays to output. Close the netcdf file and return.
2880	IF ( i == 1 ) RETURN
2881
2882	!
2883	!-- Write the list of variables as global attribute (this is used by
2884	!-- restart runs and by combine_plot_fields)
2885	nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'VAR_LIST', &
2886	var_list )
2887	CALL handle_netcdf_error( 'netcdf', 200 )
2888
2889	!
2890	!-- Leave NetCDF define mode
2891	nc_stat = NF90_ENDDEF( id_set_yz(av) )
2892	CALL handle_netcdf_error( 'netcdf', 201 )
2893
2894	!
2895	!-- Write axis data: x_yz, y, zu, zw
2896	ALLOCATE( netcdf_data(1:ns) )
2897
2898	!
2899	!-- Write x_yz data (shifted by +dx/2)
2900	DO i = 1, ns
2901	IF( section(i,3) == -1 ) THEN
2902	netcdf_data(i) = -1.0 ! section averaged along x
2903	ELSE
2904	netcdf_data(i) = ( section(i,3) + 0.5 ) * dx
2905	ENDIF
2906	ENDDO
2907	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_x_yz(av), netcdf_data, &
2908	start = (/ 1 /), count = (/ ns /) )
2909	CALL handle_netcdf_error( 'netcdf', 202 )
2910
2911	!
2912	!-- Write x_yz data (xu grid)
2913	DO i = 1, ns
2914	IF( section(i,3) == -1 ) THEN
2915	netcdf_data(i) = -1.0 ! section averaged along x
2916	ELSE
2917	netcdf_data(i) = section(i,3) * dx
2918	ENDIF
2919	ENDDO
2920	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_xu_yz(av), netcdf_data, &
2921	start = (/ 1 /), count = (/ ns /) )
2922	CALL handle_netcdf_error( 'netcdf', 383 )
2923
2924	!
2925	!-- Write gridpoint number data
2926	netcdf_data(1:ns) = section(1:ns,3)
2927	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_ind_x_yz(av), &
2928	netcdf_data, start = (/ 1 /), &
2929	count = (/ ns /) )
2930	CALL handle_netcdf_error( 'netcdf', 203 )
2931
2932	DEALLOCATE( netcdf_data )
2933
2934	!
2935	!-- Write data for y (shifted by +dy/2) and yv axis
2936	ALLOCATE( netcdf_data(0:ny+1) )
2937
2938	DO j = 0, ny+1
2939	netcdf_data(j) = ( j + 0.5 ) * dy
2940	ENDDO
2941
2942	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_y_yz(av), netcdf_data, &
2943	start = (/ 1 /), count = (/ ny+2 /) )
2944	CALL handle_netcdf_error( 'netcdf', 204 )
2945
2946	DO j = 0, ny+1
2947	netcdf_data(j) = j * dy
2948	ENDDO
2949
2950	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_yv_yz(av), &
2951	netcdf_data, start = (/ 1 /), &
2952	count = (/ ny+2 /) )
2953	CALL handle_netcdf_error( 'netcdf', 384 )
2954
2955	DEALLOCATE( netcdf_data )
2956
2957	!
2958	!-- Write zu and zw data (vertical axes)
2959	ALLOCATE( netcdf_data(0:nz+1) )
2960
2961	netcdf_data(0:nz+1) = zu(nzb:nzt+1)
2962	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zu_yz(av), &
2963	netcdf_data, start = (/ 1 /), &
2964	count = (/ nz+2 /) )
2965	CALL handle_netcdf_error( 'netcdf', 205 )
2966
2967	netcdf_data(0:nz+1) = zw(nzb:nzt+1)
2968	nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zw_yz(av), &
2969	netcdf_data, start = (/ 1 /), &
2970	count = (/ nz+2 /) )
2971	CALL handle_netcdf_error( 'netcdf', 206 )
2972
2973	DEALLOCATE( netcdf_data )
2974
2975
2976	CASE ( 'yz_ext' )
2977
2978	!
2979	!-- Get the list of variables and compare with the actual run.
2980	!-- First var_list_old has to be reset, since GET_ATT does not assign
2981	!-- trailing blanks.
2982	var_list_old = ' '
2983	nc_stat = NF90_GET_ATT( id_set_yz(av), NF90_GLOBAL, 'VAR_LIST', &
2984	var_list_old )
2985	CALL handle_netcdf_error( 'netcdf', 207 )
2986
2987	var_list = ';'
2988	i = 1
2989	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2990	IF ( INDEX( do2d(av,i), 'yz' ) /= 0 ) THEN
2991	netcdf_var_name = do2d(av,i)
2992	CALL clean_netcdf_varname( netcdf_var_name )
2993	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
2994	ENDIF
2995	i = i + 1
2996	ENDDO
2997
2998	IF ( av == 0 ) THEN
2999	var = '(yz)'
3000	ELSE
3001	var = '(yz_av)'
3002	ENDIF
3003
3004	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
3005	message_string = 'NetCDF file for cross-sections ' // &
3006	TRIM( var ) // ' from previous run found,' // &
3007	'& but this file cannot be extended due to' // &
3008	' variable mismatch.' // &
3009	'&New file is created instead.'
3010	CALL message( 'define_netcdf_header', 'PA0249', 0, 1, 0, 6, 0 )
3011	extend = .FALSE.
3012	RETURN
3013	ENDIF
3014
3015	!
3016	!-- Calculate the number of current sections
3017	ns = 1
3018	DO WHILE ( section(ns,3) /= -9999 .AND. ns <= 100 )
3019	ns = ns + 1
3020	ENDDO
3021	ns = ns - 1
3022
3023	!
3024	!-- Get and compare the number of vertical cross sections
3025	nc_stat = NF90_INQ_VARID( id_set_yz(av), 'x_yz', id_var_x_yz(av) )
3026	CALL handle_netcdf_error( 'netcdf', 208 )
3027
3028	nc_stat = NF90_INQUIRE_VARIABLE( id_set_yz(av), id_var_x_yz(av), &
3029	dimids = id_dim_x_yz_old )
3030	CALL handle_netcdf_error( 'netcdf', 209 )
3031	id_dim_x_yz(av) = id_dim_x_yz_old(1)
3032
3033	nc_stat = NF90_INQUIRE_DIMENSION( id_set_yz(av), id_dim_x_yz(av), &
3034	len = ns_old )
3035	CALL handle_netcdf_error( 'netcdf', 210 )
3036
3037	IF ( ns /= ns_old ) THEN
3038	message_string = 'NetCDF file for cross-sections ' // &
3039	TRIM( var ) // ' from previous run found,' // &
3040	'&but this file cannot be extended due to' // &
3041	' mismatch in number of' // &
3042	'&cross sections.' // &
3043	'&New file is created instead.'
3044	CALL message( 'define_netcdf_header', 'PA0250', 0, 1, 0, 6, 0 )
3045	extend = .FALSE.
3046	RETURN
3047	ENDIF
3048
3049	!
3050	!-- Get and compare the heights of the cross sections
3051	ALLOCATE( netcdf_data(1:ns_old) )
3052
3053	nc_stat = NF90_GET_VAR( id_set_yz(av), id_var_x_yz(av), netcdf_data )
3054	CALL handle_netcdf_error( 'netcdf', 211 )
3055
3056	DO i = 1, ns
3057	IF ( section(i,3) /= -1 ) THEN
3058	IF ( ( section(i,3) * dx ) /= netcdf_data(i) ) THEN
3059	message_string = 'NetCDF file for cross-sections ' // &
3060	TRIM( var ) // ' from previous run found,' // &
3061	'&but this file cannot be extended' // &
3062	' due to mismatch in cross' // &
3063	'&section levels.' // &
3064	'&New file is created instead.'
3065	CALL message( 'define_netcdf_header', 'PA0251', &
3066	0, 1, 0, 6, 0 )
3067	extend = .FALSE.
3068	RETURN
3069	ENDIF
3070	ELSE
3071	IF ( -1.0 /= netcdf_data(i) ) THEN
3072	message_string = 'NetCDF file for cross-sections ' // &
3073	TRIM( var ) // ' from previous run found,' // &
3074	'&but this file cannot be extended' // &
3075	' due to mismatch in cross' // &
3076	'&section levels.' // &
3077	'&New file is created instead.'
3078	CALL message( 'define_netcdf_header', 'PA0251', &
3079	0, 1, 0, 6, 0 )
3080	extend = .FALSE.
3081	RETURN
3082	ENDIF
3083	ENDIF
3084	ENDDO
3085
3086	DEALLOCATE( netcdf_data )
3087
3088	!
3089	!-- Get the id of the time coordinate (unlimited coordinate) and its
3090	!-- last index on the file. The next time level is pl2d..count+1.
3091	!-- The current time must be larger than the last output time
3092	!-- on the file.
3093	nc_stat = NF90_INQ_VARID( id_set_yz(av), 'time', id_var_time_yz(av) )
3094	CALL handle_netcdf_error( 'netcdf', 212 )
3095
3096	nc_stat = NF90_INQUIRE_VARIABLE( id_set_yz(av), id_var_time_yz(av), &
3097	dimids = id_dim_time_old )
3098	CALL handle_netcdf_error( 'netcdf', 213 )
3099	id_dim_time_yz(av) = id_dim_time_old(1)
3100
3101	nc_stat = NF90_INQUIRE_DIMENSION( id_set_yz(av), id_dim_time_yz(av), &
3102	len = do2d_yz_time_count(av) )
3103	CALL handle_netcdf_error( 'netcdf', 214 )
3104
3105	nc_stat = NF90_GET_VAR( id_set_yz(av), id_var_time_yz(av), &
3106	last_time_coordinate, &
3107	start = (/ do2d_yz_time_count(av) /), &
3108	count = (/ 1 /) )
3109	CALL handle_netcdf_error( 'netcdf', 215 )
3110
3111	IF ( last_time_coordinate(1) >= simulated_time ) THEN
3112	message_string = 'NetCDF file for cross sections ' // &
3113	TRIM( var ) // ' from previous run found,' // &
3114	'&but this file cannot be extended becaus' // &
3115	'e the current output time' // &
3116	'&is less or equal than the last output t' // &
3117	'ime on this file.' // &
3118	'&New file is created instead.'
3119	CALL message( 'define_netcdf_header', 'PA0252', 0, 1, 0, 6, 0 )
3120	do2d_yz_time_count(av) = 0
3121	extend = .FALSE.
3122	RETURN
3123	ENDIF
3124
3125	!
3126	!-- Dataset seems to be extendable.
3127	!-- Now get the variable ids.
3128	i = 1
3129	DO WHILE ( do2d(av,i)(1:1) /= ' ' )
3130	IF ( INDEX( do2d(av,i), 'yz' ) /= 0 ) THEN
3131	netcdf_var_name = do2d(av,i)
3132	CALL clean_netcdf_varname( netcdf_var_name )
3133	nc_stat = NF90_INQ_VARID( id_set_yz(av), netcdf_var_name, &
3134	id_var_do2d(av,i) )
3135	CALL handle_netcdf_error( 'netcdf', 216 )
3136	#if defined( __netcdf4 )
3137	!
3138	!-- Set independent io operations for parallel io. Collective io
3139	!-- is only allowed in case of a 1d-decomposition along y, because
3140	!-- otherwise, not all PEs have output data.
3141	IF ( netcdf_data_format > 2 ) THEN
3142	IF ( npex == 1 ) THEN
3143	nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av), &
3144	id_var_do2d(av,i), &
3145	NF90_COLLECTIVE )
3146	ELSE
3147	!
3148	!-- ATTENTION: Due to a probable bug in the NetCDF4
3149	!-- installation, independet output does not work
3150	!-- A workaround is used in data_output_2d on those
3151	!-- PEs having no data
3152	nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av), &
3153	id_var_do2d(av,i), &
3154	NF90_COLLECTIVE )
3155	! nc_stat = NF90_VAR_PAR_ACCESS( id_set_yz(av), &
3156	! id_var_do2d(av,i), &
3157	! NF90_INDEPENDENT )
3158	ENDIF
3159	CALL handle_netcdf_error( 'netcdf', 450 )
3160	ENDIF
3161	#endif
3162	ENDIF
3163	i = i + 1
3164	ENDDO
3165
3166	!
3167	!-- Update the title attribute on file
3168	!-- In order to avoid 'data mode' errors if updated attributes are larger
3169	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
3170	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
3171	!-- performance loss due to data copying; an alternative strategy would be
3172	!-- to ensure equal attribute size in a job chain. Maybe revise later.
3173	IF ( av == 0 ) THEN
3174	time_average_text = ' '
3175	ELSE
3176	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
3177	averaging_interval
3178	ENDIF
3179	nc_stat = NF90_REDEF( id_set_yz(av) )
3180	CALL handle_netcdf_error( 'netcdf', 435 )
3181	nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'title', &
3182	TRIM( run_description_header ) // &
3183	TRIM( time_average_text ) )
3184	CALL handle_netcdf_error( 'netcdf', 217 )
3185	nc_stat = NF90_ENDDEF( id_set_yz(av) )
3186	CALL handle_netcdf_error( 'netcdf', 436 )
3187	message_string = 'NetCDF file for cross-sections ' // &
3188	TRIM( var ) // ' from previous run found.' // &
3189	'&This file will be extended.'
3190	CALL message( 'define_netcdf_header', 'PA0253', 0, 0, 0, 6, 0 )
3191
3192
3193	CASE ( 'pr_new' )
3194
3195	!
3196	!-- Define some global attributes of the dataset
3197	IF ( averaging_interval_pr /= 0.0 ) THEN
3198	WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
3199	averaging_interval_pr
3200	nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'title', &
3201	TRIM( run_description_header ) // &
3202	TRIM( time_average_text ) )
3203	CALL handle_netcdf_error( 'netcdf', 218 )
3204
3205	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval_pr
3206	nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'time_avg', &
3207	TRIM( time_average_text ) )
3208	ELSE
3209	nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'title', &
3210	TRIM( run_description_header ) )
3211	ENDIF
3212	CALL handle_netcdf_error( 'netcdf', 219 )
3213
3214	!
3215	!-- Define time coordinate for profiles (unlimited dimension)
3216	nc_stat = NF90_DEF_DIM( id_set_pr, 'time', NF90_UNLIMITED, &
3217	id_dim_time_pr )
3218	CALL handle_netcdf_error( 'netcdf', 220 )
3219
3220	nc_stat = NF90_DEF_VAR( id_set_pr, 'time', NF90_DOUBLE, &
3221	id_dim_time_pr, id_var_time_pr )
3222	CALL handle_netcdf_error( 'netcdf', 221 )
3223
3224	nc_stat = NF90_PUT_ATT( id_set_pr, id_var_time_pr, 'units', 'seconds')
3225	CALL handle_netcdf_error( 'netcdf', 222 )
3226
3227	!
3228	!-- Define the variables
3229	var_list = ';'
3230	DO i = 1, dopr_n
3231	!
3232	!-- First, remove those characters not allowed by NetCDF
3233	netcdf_var_name = data_output_pr(i)
3234	CALL clean_netcdf_varname( netcdf_var_name )
3235
3236	IF ( statistic_regions == 0 ) THEN
3237
3238	!
3239	!-- Define the z-axes (each variable gets its own z-axis)
3240	nc_stat = NF90_DEF_DIM( id_set_pr, 'z'//TRIM(netcdf_var_name), &
3241	nzt+2-nzb, id_dim_z_pr(i,0) )
3242	CALL handle_netcdf_error( 'netcdf', 223 )
3243
3244	nc_stat = NF90_DEF_VAR( id_set_pr, 'z'//TRIM(netcdf_var_name), &
3245	NF90_DOUBLE, id_dim_z_pr(i,0), &
3246	id_var_z_pr(i,0) )
3247	CALL handle_netcdf_error( 'netcdf', 224 )
3248
3249	nc_stat = NF90_PUT_ATT( id_set_pr, id_var_z_pr(i,0), 'units', &
3250	'meters' )
3251	CALL handle_netcdf_error( 'netcdf', 225 )
3252	!
3253	!-- Define the variable
3254	nc_stat = NF90_DEF_VAR( id_set_pr, netcdf_var_name, &
3255	nc_precision(5), (/ id_dim_z_pr(i,0), &
3256	id_dim_time_pr /), id_var_dopr(i,0) )
3257	CALL handle_netcdf_error( 'netcdf', 226 )
3258
3259	nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,0), &
3260	'long_name', TRIM( data_output_pr(i) ) )
3261	CALL handle_netcdf_error( 'netcdf', 227 )
3262
3263	nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,0), &
3264	'units', TRIM( dopr_unit(i) ) )
3265	CALL handle_netcdf_error( 'netcdf', 228 )
3266
3267	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
3268
3269	ELSE
3270	!
3271	!-- If statistic regions are defined, add suffix _SR+#SR to the
3272	!-- names
3273	DO j = 0, statistic_regions
3274	WRITE ( suffix, '(''_'',I1)' ) j
3275
3276	!
3277	!-- Define the z-axes (each variable gets it own z-axis)
3278	nc_stat = NF90_DEF_DIM( id_set_pr, &
3279	'z'//TRIM(netcdf_var_name)//suffix, &
3280	nzt+2-nzb, id_dim_z_pr(i,j) )
3281	CALL handle_netcdf_error( 'netcdf', 229 )
3282
3283	nc_stat = NF90_DEF_VAR( id_set_pr, &
3284	'z'//TRIM(netcdf_var_name)//suffix, &
3285	nc_precision(5), id_dim_z_pr(i,j), &
3286	id_var_z_pr(i,j) )
3287	CALL handle_netcdf_error( 'netcdf', 230 )
3288
3289	nc_stat = NF90_PUT_ATT( id_set_pr, id_var_z_pr(i,j), &
3290	'units', 'meters' )
3291	CALL handle_netcdf_error( 'netcdf', 231 )
3292
3293	!
3294	!-- Define the variable
3295	nc_stat = NF90_DEF_VAR( id_set_pr, &
3296	TRIM( netcdf_var_name ) // suffix, &
3297	nc_precision(5), &
3298	(/ id_dim_z_pr(i,j), &
3299	id_dim_time_pr /), id_var_dopr(i,j) )
3300	CALL handle_netcdf_error( 'netcdf', 232 )
3301
3302	nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,j), &
3303	'long_name', &
3304	TRIM( data_output_pr(i) ) // ' SR ' &
3305	// suffix(2:2) )
3306	CALL handle_netcdf_error( 'netcdf', 233 )
3307
3308	nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,j), &
3309	'units', TRIM( dopr_unit(i) ) )
3310	CALL handle_netcdf_error( 'netcdf', 234 )
3311
3312	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
3313	suffix // ';'
3314
3315	ENDDO
3316
3317	ENDIF
3318
3319	ENDDO
3320
3321	!
3322	!-- Write the list of variables as global attribute (this is used by
3323	!-- restart runs)
3324	nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'VAR_LIST', var_list )
3325	CALL handle_netcdf_error( 'netcdf', 235 )
3326
3327	!
3328	!-- Define normalization variables (as time series)
3329	DO i = 1, dopr_norm_num
3330
3331	nc_stat = NF90_DEF_VAR( id_set_pr, 'NORM_' // &
3332	TRIM( dopr_norm_names(i) ), &
3333	nc_precision(5), (/ id_dim_time_pr /), &
3334	id_var_norm_dopr(i) )
3335	CALL handle_netcdf_error( 'netcdf', 236 )
3336
3337	nc_stat = NF90_PUT_ATT( id_set_pr, id_var_norm_dopr(i), &
3338	'long_name', &
3339	TRIM( dopr_norm_longnames(i) ) )
3340	CALL handle_netcdf_error( 'netcdf', 237 )
3341
3342	ENDDO
3343
3344	!
3345	!-- Leave NetCDF define mode
3346	nc_stat = NF90_ENDDEF( id_set_pr )
3347	CALL handle_netcdf_error( 'netcdf', 238 )
3348
3349	!
3350	!-- Write z-axes data
3351	DO i = 1, dopr_n
3352	DO j = 0, statistic_regions
3353
3354	nc_stat = NF90_PUT_VAR( id_set_pr, id_var_z_pr(i,j), &
3355	hom(nzb:nzt+1,2,dopr_index(i),0), &
3356	start = (/ 1 /), &
3357	count = (/ nzt-nzb+2 /) )
3358	CALL handle_netcdf_error( 'netcdf', 239 )
3359
3360	ENDDO
3361	ENDDO
3362
3363
3364	CASE ( 'pr_ext' )
3365
3366	!
3367	!-- Get the list of variables and compare with the actual run.
3368	!-- First var_list_old has to be reset, since GET_ATT does not assign
3369	!-- trailing blanks.
3370	var_list_old = ' '
3371	nc_stat = NF90_GET_ATT( id_set_pr, NF90_GLOBAL, 'VAR_LIST', &
3372	var_list_old )
3373	CALL handle_netcdf_error( 'netcdf', 240 )
3374
3375	var_list = ';'
3376	DO i = 1, dopr_n
3377
3378	netcdf_var_name = data_output_pr(i)
3379	CALL clean_netcdf_varname( netcdf_var_name )
3380
3381	IF ( statistic_regions == 0 ) THEN
3382	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
3383	ELSE
3384	DO j = 0, statistic_regions
3385	WRITE ( suffix, '(''_'',I1)' ) j
3386	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
3387	suffix // ';'
3388	ENDDO
3389	ENDIF
3390
3391	ENDDO
3392
3393	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
3394	message_string = 'NetCDF file for vertical profiles ' // &
3395	'from previous run found,' // &
3396	'& but this file cannot be extended due to' // &
3397	' variable mismatch.' // &
3398	'&New file is created instead.'
3399	CALL message( 'define_netcdf_header', 'PA0254', 0, 1, 0, 6, 0 )
3400	extend = .FALSE.
3401	RETURN
3402	ENDIF
3403
3404	!
3405	!-- Get the id of the time coordinate (unlimited coordinate) and its
3406	!-- last index on the file. The next time level is dopr..count+1.
3407	!-- The current time must be larger than the last output time
3408	!-- on the file.
3409	nc_stat = NF90_INQ_VARID( id_set_pr, 'time', id_var_time_pr )
3410	CALL handle_netcdf_error( 'netcdf', 241 )
3411
3412	nc_stat = NF90_INQUIRE_VARIABLE( id_set_pr, id_var_time_pr, &
3413	dimids = id_dim_time_old )
3414	CALL handle_netcdf_error( 'netcdf', 242 )
3415	id_dim_time_pr = id_dim_time_old(1)
3416
3417	nc_stat = NF90_INQUIRE_DIMENSION( id_set_pr, id_dim_time_pr, &
3418	len = dopr_time_count )
3419	CALL handle_netcdf_error( 'netcdf', 243 )
3420
3421	nc_stat = NF90_GET_VAR( id_set_pr, id_var_time_pr, &
3422	last_time_coordinate, &
3423	start = (/ dopr_time_count /), &
3424	count = (/ 1 /) )
3425	CALL handle_netcdf_error( 'netcdf', 244 )
3426
3427	IF ( last_time_coordinate(1) >= simulated_time ) THEN
3428	message_string = 'NetCDF file for vertical profiles ' // &
3429	'from previous run found,' // &
3430	'&but this file cannot be extended becaus' // &
3431	'e the current output time' // &
3432	'&is less or equal than the last output t' // &
3433	'ime on this file.' // &
3434	'&New file is created instead.'
3435	CALL message( 'define_netcdf_header', 'PA0255', 0, 1, 0, 6, 0 )
3436	dopr_time_count = 0
3437	extend = .FALSE.
3438	RETURN
3439	ENDIF
3440
3441	!
3442	!-- Dataset seems to be extendable.
3443	!-- Now get the variable ids.
3444	i = 1
3445	DO i = 1, dopr_n
3446
3447	netcdf_var_name_base = data_output_pr(i)
3448	CALL clean_netcdf_varname( netcdf_var_name_base )
3449
3450	IF ( statistic_regions == 0 ) THEN
3451	nc_stat = NF90_INQ_VARID( id_set_pr, netcdf_var_name_base, &
3452	id_var_dopr(i,0) )
3453	CALL handle_netcdf_error( 'netcdf', 245 )
3454	ELSE
3455	DO j = 0, statistic_regions
3456	WRITE ( suffix, '(''_'',I1)' ) j
3457	netcdf_var_name = TRIM( netcdf_var_name_base ) // suffix
3458	nc_stat = NF90_INQ_VARID( id_set_pr, netcdf_var_name, &
3459	id_var_dopr(i,j) )
3460	CALL handle_netcdf_error( 'netcdf', 246 )
3461	ENDDO
3462	ENDIF
3463
3464	ENDDO
3465
3466	!
3467	!-- Get ids of the normalization variables
3468	DO i = 1, dopr_norm_num
3469	nc_stat = NF90_INQ_VARID( id_set_pr, &
3470	'NORM_' // TRIM( dopr_norm_names(i) ), &
3471	id_var_norm_dopr(i) )
3472	CALL handle_netcdf_error( 'netcdf', 247 )
3473	ENDDO
3474
3475	!
3476	!-- Update the title attribute on file
3477	!-- In order to avoid 'data mode' errors if updated attributes are larger
3478	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
3479	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
3480	!-- performance loss due to data copying; an alternative strategy would be
3481	!-- to ensure equal attribute size in a job chain. Maybe revise later.
3482	IF ( averaging_interval_pr == 0.0 ) THEN
3483	time_average_text = ' '
3484	ELSE
3485	WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
3486	averaging_interval_pr
3487	ENDIF
3488	nc_stat = NF90_REDEF( id_set_pr )
3489	CALL handle_netcdf_error( 'netcdf', 437 )
3490	nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'title', &
3491	TRIM( run_description_header ) // &
3492	TRIM( time_average_text ) )
3493	CALL handle_netcdf_error( 'netcdf', 248 )
3494	nc_stat = NF90_ENDDEF( id_set_pr )
3495	CALL handle_netcdf_error( 'netcdf', 438 )
3496	message_string = 'NetCDF file for vertical profiles ' // &
3497	'from previous run found.' // &
3498	'&This file will be extended.'
3499	CALL message( 'define_netcdf_header', 'PA0256', 0, 0, 0, 6, 0 )
3500
3501
3502	CASE ( 'ts_new' )
3503
3504	!
3505	!-- Define some global attributes of the dataset
3506	nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'title', &
3507	TRIM( run_description_header ) )
3508	CALL handle_netcdf_error( 'netcdf', 249 )
3509
3510	!
3511	!-- Define time coordinate for time series (unlimited dimension)
3512	nc_stat = NF90_DEF_DIM( id_set_ts, 'time', NF90_UNLIMITED, &
3513	id_dim_time_ts )
3514	CALL handle_netcdf_error( 'netcdf', 250 )
3515
3516	nc_stat = NF90_DEF_VAR( id_set_ts, 'time', NF90_DOUBLE, &
3517	id_dim_time_ts, id_var_time_ts )
3518	CALL handle_netcdf_error( 'netcdf', 251 )
3519
3520	nc_stat = NF90_PUT_ATT( id_set_ts, id_var_time_ts, 'units', 'seconds')
3521	CALL handle_netcdf_error( 'netcdf', 252 )
3522
3523	!
3524	!-- Define the variables
3525	var_list = ';'
3526	DO i = 1, dots_num
3527	!
3528	!-- First, remove those characters not allowed by NetCDF
3529	netcdf_var_name = dots_label(i)
3530	CALL clean_netcdf_varname( netcdf_var_name )
3531
3532	IF ( statistic_regions == 0 ) THEN
3533
3534	nc_stat = NF90_DEF_VAR( id_set_ts, netcdf_var_name, &
3535	nc_precision(6), (/ id_dim_time_ts /), &
3536	id_var_dots(i,0) )
3537	CALL handle_netcdf_error( 'netcdf', 253 )
3538
3539	nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,0), &
3540	'long_name', TRIM( dots_label(i) ) )
3541	CALL handle_netcdf_error( 'netcdf', 254 )
3542
3543	nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,0), &
3544	'units', TRIM( dots_unit(i) ) )
3545	CALL handle_netcdf_error( 'netcdf', 255 )
3546
3547	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
3548
3549	ELSE
3550	!
3551	!-- If statistic regions are defined, add suffix _SR+#SR to the
3552	!-- names
3553	DO j = 0, statistic_regions
3554	WRITE ( suffix, '(''_'',I1)' ) j
3555
3556	nc_stat = NF90_DEF_VAR( id_set_ts, &
3557	TRIM( netcdf_var_name ) // suffix, &
3558	nc_precision(6), &
3559	(/ id_dim_time_ts /), &
3560	id_var_dots(i,j) )
3561	CALL handle_netcdf_error( 'netcdf', 256 )
3562
3563	nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,j), &
3564	'long_name', &
3565	TRIM( dots_label(i) ) // ' SR ' // &
3566	suffix(2:2) )
3567	CALL handle_netcdf_error( 'netcdf', 257 )
3568
3569	nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,j), &
3570	'units', TRIM( dots_unit(i) ) )
3571	CALL handle_netcdf_error( 'netcdf', 347 )
3572
3573	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
3574	suffix // ';'
3575
3576	ENDDO
3577
3578	ENDIF
3579
3580	ENDDO
3581
3582	!
3583	!-- Write the list of variables as global attribute (this is used by
3584	!-- restart runs)
3585	nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'VAR_LIST', var_list )
3586	CALL handle_netcdf_error( 'netcdf', 258 )
3587
3588	!
3589	!-- Leave NetCDF define mode
3590	nc_stat = NF90_ENDDEF( id_set_ts )
3591	CALL handle_netcdf_error( 'netcdf', 259 )
3592
3593
3594	CASE ( 'ts_ext' )
3595
3596	!
3597	!-- Get the list of variables and compare with the actual run.
3598	!-- First var_list_old has to be reset, since GET_ATT does not assign
3599	!-- trailing blanks.
3600	var_list_old = ' '
3601	nc_stat = NF90_GET_ATT( id_set_ts, NF90_GLOBAL, 'VAR_LIST', &
3602	var_list_old )
3603	CALL handle_netcdf_error( 'netcdf', 260 )
3604
3605	var_list = ';'
3606	i = 1
3607	DO i = 1, dots_num
3608
3609	netcdf_var_name = dots_label(i)
3610	CALL clean_netcdf_varname( netcdf_var_name )
3611
3612	IF ( statistic_regions == 0 ) THEN
3613	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
3614	ELSE
3615	DO j = 0, statistic_regions
3616	WRITE ( suffix, '(''_'',I1)' ) j
3617	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
3618	suffix // ';'
3619	ENDDO
3620	ENDIF
3621
3622	ENDDO
3623
3624	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
3625	message_string = 'NetCDF file for time series ' // &
3626	'from previous run found,' // &
3627	'& but this file cannot be extended due to' // &
3628	' variable mismatch.' // &
3629	'&New file is created instead.'
3630	CALL message( 'define_netcdf_header', 'PA0257', 0, 1, 0, 6, 0 )
3631	extend = .FALSE.
3632	RETURN
3633	ENDIF
3634
3635	!
3636	!-- Get the id of the time coordinate (unlimited coordinate) and its
3637	!-- last index on the file. The next time level is dots..count+1.
3638	!-- The current time must be larger than the last output time
3639	!-- on the file.
3640	nc_stat = NF90_INQ_VARID( id_set_ts, 'time', id_var_time_ts )
3641	CALL handle_netcdf_error( 'netcdf', 261 )
3642
3643	nc_stat = NF90_INQUIRE_VARIABLE( id_set_ts, id_var_time_ts, &
3644	dimids = id_dim_time_old )
3645	CALL handle_netcdf_error( 'netcdf', 262 )
3646	id_dim_time_ts = id_dim_time_old(1)
3647
3648	nc_stat = NF90_INQUIRE_DIMENSION( id_set_ts, id_dim_time_ts, &
3649	len = dots_time_count )
3650	CALL handle_netcdf_error( 'netcdf', 263 )
3651
3652	nc_stat = NF90_GET_VAR( id_set_ts, id_var_time_ts, &
3653	last_time_coordinate, &
3654	start = (/ dots_time_count /), &
3655	count = (/ 1 /) )
3656	CALL handle_netcdf_error( 'netcdf', 264 )
3657
3658	IF ( last_time_coordinate(1) >= simulated_time ) THEN
3659	message_string = 'NetCDF file for time series ' // &
3660	'from previous run found,' // &
3661	'&but this file cannot be extended becaus' // &
3662	'e the current output time' // &
3663	'&is less or equal than the last output t' // &
3664	'ime on this file.' // &
3665	'&New file is created instead.'
3666	CALL message( 'define_netcdf_header', 'PA0258', 0, 1, 0, 6, 0 )
3667	dots_time_count = 0
3668	extend = .FALSE.
3669	RETURN
3670	ENDIF
3671
3672	!
3673	!-- Dataset seems to be extendable.
3674	!-- Now get the variable ids
3675	i = 1
3676	DO i = 1, dots_num
3677
3678	netcdf_var_name_base = dots_label(i)
3679	CALL clean_netcdf_varname( netcdf_var_name_base )
3680
3681	IF ( statistic_regions == 0 ) THEN
3682	nc_stat = NF90_INQ_VARID( id_set_ts, netcdf_var_name_base, &
3683	id_var_dots(i,0) )
3684	CALL handle_netcdf_error( 'netcdf', 265 )
3685	ELSE
3686	DO j = 0, statistic_regions
3687	WRITE ( suffix, '(''_'',I1)' ) j
3688	netcdf_var_name = TRIM( netcdf_var_name_base ) // suffix
3689	nc_stat = NF90_INQ_VARID( id_set_ts, netcdf_var_name, &
3690	id_var_dots(i,j) )
3691	CALL handle_netcdf_error( 'netcdf', 266 )
3692	ENDDO
3693	ENDIF
3694
3695	ENDDO
3696
3697	!
3698	!-- Update the title attribute on file
3699	!-- In order to avoid 'data mode' errors if updated attributes are larger
3700	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
3701	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
3702	!-- performance loss due to data copying; an alternative strategy would be
3703	!-- to ensure equal attribute size in a job chain. Maybe revise later.
3704	nc_stat = NF90_REDEF( id_set_ts )
3705	CALL handle_netcdf_error( 'netcdf', 439 )
3706	nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'title', &
3707	TRIM( run_description_header ) )
3708	CALL handle_netcdf_error( 'netcdf', 267 )
3709	nc_stat = NF90_ENDDEF( id_set_ts )
3710	CALL handle_netcdf_error( 'netcdf', 440 )
3711	message_string = 'NetCDF file for time series ' // &
3712	'from previous run found.' // &
3713	'&This file will be extended.'
3714	CALL message( 'define_netcdf_header', 'PA0259', 0, 0, 0, 6, 0 )
3715
3716
3717	CASE ( 'sp_new' )
3718
3719	!
3720	!-- Define some global attributes of the dataset
3721	IF ( averaging_interval_sp /= 0.0 ) THEN
3722	WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
3723	averaging_interval_sp
3724	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
3725	TRIM( run_description_header ) // &
3726	TRIM( time_average_text ) )
3727	CALL handle_netcdf_error( 'netcdf', 268 )
3728
3729	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval_sp
3730	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'time_avg', &
3731	TRIM( time_average_text ) )
3732	ELSE
3733	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
3734	TRIM( run_description_header ) )
3735	ENDIF
3736	CALL handle_netcdf_error( 'netcdf', 269 )
3737
3738	!
3739	!-- Define time coordinate for spectra (unlimited dimension)
3740	nc_stat = NF90_DEF_DIM( id_set_sp, 'time', NF90_UNLIMITED, &
3741	id_dim_time_sp )
3742	CALL handle_netcdf_error( 'netcdf', 270 )
3743
3744	nc_stat = NF90_DEF_VAR( id_set_sp, 'time', NF90_DOUBLE, &
3745	id_dim_time_sp, id_var_time_sp )
3746	CALL handle_netcdf_error( 'netcdf', 271 )
3747
3748	nc_stat = NF90_PUT_ATT( id_set_sp, id_var_time_sp, 'units', 'seconds')
3749	CALL handle_netcdf_error( 'netcdf', 272 )
3750
3751	!
3752	!-- Define the spatial dimensions and coordinates for spectra.
3753	!-- First, determine the number of vertical levels for which spectra
3754	!-- are to be output.
3755	ns = 1
3756	DO WHILE ( comp_spectra_level(ns) /= 999999 .AND. ns <= 100 )
3757	ns = ns + 1
3758	ENDDO
3759	ns = ns - 1
3760
3761	!
3762	!-- Define vertical coordinate grid (zu grid)
3763	nc_stat = NF90_DEF_DIM( id_set_sp, 'zu_sp', ns, id_dim_zu_sp )
3764	CALL handle_netcdf_error( 'netcdf', 273 )
3765
3766	nc_stat = NF90_DEF_VAR( id_set_sp, 'zu_sp', NF90_DOUBLE, &
3767	id_dim_zu_sp, id_var_zu_sp )
3768	CALL handle_netcdf_error( 'netcdf', 274 )
3769
3770	nc_stat = NF90_PUT_ATT( id_set_sp, id_var_zu_sp, 'units', 'meters' )
3771	CALL handle_netcdf_error( 'netcdf', 275 )
3772
3773	!
3774	!-- Define vertical coordinate grid (zw grid)
3775	nc_stat = NF90_DEF_DIM( id_set_sp, 'zw_sp', ns, id_dim_zw_sp )
3776	CALL handle_netcdf_error( 'netcdf', 276 )
3777
3778	nc_stat = NF90_DEF_VAR( id_set_sp, 'zw_sp', NF90_DOUBLE, &
3779	id_dim_zw_sp, id_var_zw_sp )
3780	CALL handle_netcdf_error( 'netcdf', 277 )
3781
3782	nc_stat = NF90_PUT_ATT( id_set_sp, id_var_zw_sp, 'units', 'meters' )
3783	CALL handle_netcdf_error( 'netcdf', 278 )
3784
3785	!
3786	!-- Define x-axis
3787	nc_stat = NF90_DEF_DIM( id_set_sp, 'k_x', nx/2, id_dim_x_sp )
3788	CALL handle_netcdf_error( 'netcdf', 279 )
3789
3790	nc_stat = NF90_DEF_VAR( id_set_sp, 'k_x', NF90_DOUBLE, id_dim_x_sp, &
3791	id_var_x_sp )
3792	CALL handle_netcdf_error( 'netcdf', 280 )
3793
3794	nc_stat = NF90_PUT_ATT( id_set_sp, id_var_x_sp, 'units', 'm-1' )
3795	CALL handle_netcdf_error( 'netcdf', 281 )
3796
3797	!
3798	!-- Define y-axis
3799	nc_stat = NF90_DEF_DIM( id_set_sp, 'k_y', ny/2, id_dim_y_sp )
3800	CALL handle_netcdf_error( 'netcdf', 282 )
3801
3802	nc_stat = NF90_DEF_VAR( id_set_sp, 'k_y', NF90_DOUBLE, id_dim_y_sp, &
3803	id_var_y_sp )
3804	CALL handle_netcdf_error( 'netcdf', 283 )
3805
3806	nc_stat = NF90_PUT_ATT( id_set_sp, id_var_y_sp, 'units', 'm-1' )
3807	CALL handle_netcdf_error( 'netcdf', 284 )
3808
3809	!
3810	!-- Define the variables
3811	var_list = ';'
3812	i = 1
3813	DO WHILE ( data_output_sp(i) /= ' ' .AND. i <= 10 )
3814	!
3815	!-- First check for the vertical grid
3816	found = .TRUE.
3817	SELECT CASE ( data_output_sp(i) )
3818	!
3819	!-- Most variables are defined on the zu levels
3820	CASE ( 'e', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c', 'ql_v', &
3821	'ql_vp', 'qv', 'rho', 's', 'sa', 'u', 'v', 'vpt' )
3822
3823	grid_z = 'zu'
3824	!
3825	!-- zw levels
3826	CASE ( 'w' )
3827
3828	grid_z = 'zw'
3829
3830	CASE DEFAULT
3831	!
3832	!-- Check for user-defined quantities (found, grid_x and grid_y
3833	!-- are dummies)
3834	CALL user_define_netcdf_grid( data_output_sp(i), found, &
3835	grid_x, grid_y, grid_z )
3836
3837	END SELECT
3838
3839	IF ( INDEX( spectra_direction(i), 'x' ) /= 0 ) THEN
3840
3841	!
3842	!-- Define the variable
3843	netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
3844	CALL clean_netcdf_varname( netcdf_var_name )
3845	IF ( TRIM( grid_z ) == 'zw' ) THEN
3846	nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name, &
3847	nc_precision(7), (/ id_dim_x_sp, &
3848	id_dim_zw_sp, id_dim_time_sp /), &
3849	id_var_dospx(i) )
3850	ELSE
3851	nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name, &
3852	nc_precision(7), (/ id_dim_x_sp, &
3853	id_dim_zu_sp, id_dim_time_sp /), &
3854	id_var_dospx(i) )
3855	ENDIF
3856	CALL handle_netcdf_error( 'netcdf', 285 )
3857
3858	nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospx(i), &
3859	'long_name', netcdf_var_name )
3860	CALL handle_netcdf_error( 'netcdf', 286 )
3861
3862	nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospx(i), &
3863	'units', 'unknown' )
3864	CALL handle_netcdf_error( 'netcdf', 287 )
3865
3866	var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
3867
3868	ENDIF
3869
3870	IF ( INDEX( spectra_direction(i), 'y' ) /= 0 ) THEN
3871
3872	!
3873	!-- Define the variable
3874	netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
3875	CALL clean_netcdf_varname( netcdf_var_name )
3876	IF ( TRIM( grid_z ) == 'zw' ) THEN
3877	nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name, &
3878	nc_precision(7), (/ id_dim_y_sp, &
3879	id_dim_zw_sp, id_dim_time_sp /), &
3880	id_var_dospy(i) )
3881	ELSE
3882	nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name, &
3883	nc_precision(7), (/ id_dim_y_sp, &
3884	id_dim_zu_sp, id_dim_time_sp /), &
3885	id_var_dospy(i) )
3886	ENDIF
3887	CALL handle_netcdf_error( 'netcdf', 288 )
3888
3889	nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospy(i), &
3890	'long_name', netcdf_var_name )
3891	CALL handle_netcdf_error( 'netcdf', 289 )
3892
3893	nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospy(i), &
3894	'units', 'unknown' )
3895	CALL handle_netcdf_error( 'netcdf', 290 )
3896
3897	var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
3898
3899	ENDIF
3900
3901	i = i + 1
3902
3903	ENDDO
3904
3905	!
3906	!-- Write the list of variables as global attribute (this is used by
3907	!-- restart runs)
3908	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', var_list )
3909	CALL handle_netcdf_error( 'netcdf', 291 )
3910
3911	!
3912	!-- Leave NetCDF define mode
3913	nc_stat = NF90_ENDDEF( id_set_sp )
3914	CALL handle_netcdf_error( 'netcdf', 292 )
3915
3916	!
3917	!-- Write axis data: zu_sp, zw_sp, k_x, k_y
3918	ALLOCATE( netcdf_data(1:ns) )
3919
3920	!
3921	!-- Write zu data
3922	netcdf_data(1:ns) = zu( comp_spectra_level(1:ns) )
3923	nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zu_sp, netcdf_data, &
3924	start = (/ 1 /), count = (/ ns /) )
3925	CALL handle_netcdf_error( 'netcdf', 293 )
3926
3927	!
3928	!-- Write zw data
3929	netcdf_data(1:ns) = zw( comp_spectra_level(1:ns) )
3930	nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zw_sp, netcdf_data, &
3931	start = (/ 1 /), count = (/ ns /) )
3932	CALL handle_netcdf_error( 'netcdf', 294 )
3933
3934	DEALLOCATE( netcdf_data )
3935
3936	!
3937	!-- Write data for x and y axis (wavenumbers)
3938	ALLOCATE( netcdf_data(nx/2) )
3939	DO i = 1, nx/2
3940	netcdf_data(i) = 2.0 * pi * i / ( dx * ( nx + 1 ) )
3941	ENDDO
3942
3943	nc_stat = NF90_PUT_VAR( id_set_sp, id_var_x_sp, netcdf_data, &
3944	start = (/ 1 /), count = (/ nx/2 /) )
3945	CALL handle_netcdf_error( 'netcdf', 295 )
3946
3947	DEALLOCATE( netcdf_data )
3948
3949	ALLOCATE( netcdf_data(ny/2) )
3950	DO i = 1, ny/2
3951	netcdf_data(i) = 2.0 * pi * i / ( dy * ( ny + 1 ) )
3952	ENDDO
3953
3954	nc_stat = NF90_PUT_VAR( id_set_sp, id_var_y_sp, netcdf_data, &
3955	start = (/ 1 /), count = (/ ny/2 /) )
3956	CALL handle_netcdf_error( 'netcdf', 296 )
3957
3958	DEALLOCATE( netcdf_data )
3959
3960
3961	CASE ( 'sp_ext' )
3962
3963	!
3964	!-- Get the list of variables and compare with the actual run.
3965	!-- First var_list_old has to be reset, since GET_ATT does not assign
3966	!-- trailing blanks.
3967	var_list_old = ' '
3968	nc_stat = NF90_GET_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', &
3969	var_list_old )
3970	CALL handle_netcdf_error( 'netcdf', 297 )
3971	var_list = ';'
3972	i = 1
3973	DO WHILE ( data_output_sp(i) /= ' ' .AND. i <= 10 )
3974
3975	IF ( INDEX( spectra_direction(i), 'x' ) /= 0 ) THEN
3976	netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
3977	var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
3978	ENDIF
3979
3980	IF ( INDEX( spectra_direction(i), 'y' ) /= 0 ) THEN
3981	netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
3982	var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
3983	ENDIF
3984
3985	i = i + 1
3986
3987	ENDDO
3988
3989	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
3990	message_string = 'NetCDF file for spectra ' // &
3991	'from previous run found,' // &
3992	'& but this file cannot be extended due to' // &
3993	' variable mismatch.' // &
3994	'&New file is created instead.'
3995	CALL message( 'define_netcdf_header', 'PA0260', 0, 1, 0, 6, 0 )
3996	extend = .FALSE.
3997	RETURN
3998	ENDIF
3999
4000	!
4001	!-- Determine the number of current vertical levels for which spectra
4002	!-- shall be output
4003	ns = 1
4004	DO WHILE ( comp_spectra_level(ns) /= 999999 .AND. ns <= 100 )
4005	ns = ns + 1
4006	ENDDO
4007	ns = ns - 1
4008
4009	!
4010	!-- Get and compare the number of vertical levels
4011	nc_stat = NF90_INQ_VARID( id_set_sp, 'zu_sp', id_var_zu_sp )
4012	CALL handle_netcdf_error( 'netcdf', 298 )
4013
4014	nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_zu_sp, &
4015	dimids = id_dim_zu_sp_old )
4016	CALL handle_netcdf_error( 'netcdf', 299 )
4017	id_dim_zu_sp = id_dim_zu_sp_old(1)
4018
4019	nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_zu_sp, &
4020	len = ns_old )
4021	CALL handle_netcdf_error( 'netcdf', 300 )
4022
4023	IF ( ns /= ns_old ) THEN
4024	message_string = 'NetCDF file for spectra ' // &
4025	' from previous run found,' // &
4026	'&but this file cannot be extended due to' // &
4027	' mismatch in number of' // &
4028	'&vertical levels.' // &
4029	'&New file is created instead.'
4030	CALL message( 'define_netcdf_header', 'PA0261', 0, 1, 0, 6, 0 )
4031	extend = .FALSE.
4032	RETURN
4033	ENDIF
4034
4035	!
4036	!-- Get and compare the heights of the cross sections
4037	ALLOCATE( netcdf_data(1:ns_old) )
4038
4039	nc_stat = NF90_GET_VAR( id_set_sp, id_var_zu_sp, netcdf_data )
4040	CALL handle_netcdf_error( 'netcdf', 301 )
4041
4042	DO i = 1, ns
4043	IF ( zu(comp_spectra_level(i)) /= netcdf_data(i) ) THEN
4044	message_string = 'NetCDF file for spectra ' // &
4045	' from previous run found,' // &
4046	'&but this file cannot be extended due to' // &
4047	' mismatch in heights of' // &
4048	'&vertical levels.' // &
4049	'&New file is created instead.'
4050	CALL message( 'define_netcdf_header', 'PA0262', 0, 1, 0, 6, 0 )
4051	extend = .FALSE.
4052	RETURN
4053	ENDIF
4054	ENDDO
4055
4056	DEALLOCATE( netcdf_data )
4057
4058	!
4059	!-- Get the id of the time coordinate (unlimited coordinate) and its
4060	!-- last index on the file. The next time level is plsp..count+1.
4061	!-- The current time must be larger than the last output time
4062	!-- on the file.
4063	nc_stat = NF90_INQ_VARID( id_set_sp, 'time', id_var_time_sp )
4064	CALL handle_netcdf_error( 'netcdf', 302 )
4065
4066	nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_time_sp, &
4067	dimids = id_dim_time_old )
4068	CALL handle_netcdf_error( 'netcdf', 303 )
4069	id_dim_time_sp = id_dim_time_old(1)
4070
4071	nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_time_sp, &
4072	len = dosp_time_count )
4073	CALL handle_netcdf_error( 'netcdf', 304 )
4074
4075	nc_stat = NF90_GET_VAR( id_set_sp, id_var_time_sp, &
4076	last_time_coordinate, &
4077	start = (/ dosp_time_count /), &
4078	count = (/ 1 /) )
4079	CALL handle_netcdf_error( 'netcdf', 305 )
4080
4081	IF ( last_time_coordinate(1) >= simulated_time ) THEN
4082	message_string = 'NetCDF file for spectra ' // &
4083	'from previous run found,' // &
4084	'&but this file cannot be extended becaus' // &
4085	'e the current output time' // &
4086	'&is less or equal than the last output t' // &
4087	'ime on this file.' // &
4088	'&New file is created instead.'
4089	CALL message( 'define_netcdf_header', 'PA0263', 0, 1, 0, 6, 0 )
4090	dosp_time_count = 0
4091	extend = .FALSE.
4092	RETURN
4093	ENDIF
4094
4095	!
4096	!-- Dataset seems to be extendable.
4097	!-- Now get the variable ids.
4098	i = 1
4099	DO WHILE ( data_output_sp(i) /= ' ' .AND. i <= 10 )
4100
4101	IF ( INDEX( spectra_direction(i), 'x' ) /= 0 ) THEN
4102	netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
4103	nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
4104	id_var_dospx(i) )
4105	CALL handle_netcdf_error( 'netcdf', 306 )
4106	ENDIF
4107
4108	IF ( INDEX( spectra_direction(i), 'y' ) /= 0 ) THEN
4109	netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
4110	nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
4111	id_var_dospy(i) )
4112	CALL handle_netcdf_error( 'netcdf', 307 )
4113	ENDIF
4114
4115	i = i + 1
4116
4117	ENDDO
4118
4119	!
4120	!-- Update the title attribute on file
4121	!-- In order to avoid 'data mode' errors if updated attributes are larger
4122	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
4123	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
4124	!-- performance loss due to data copying; an alternative strategy would be
4125	!-- to ensure equal attribute size in a job chain. Maybe revise later.
4126	nc_stat = NF90_REDEF( id_set_sp )
4127	CALL handle_netcdf_error( 'netcdf', 441 )
4128	IF ( averaging_interval_sp /= 0.0 ) THEN
4129	WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
4130	averaging_interval_sp
4131	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
4132	TRIM( run_description_header ) // &
4133	TRIM( time_average_text ) )
4134	CALL handle_netcdf_error( 'netcdf', 308 )
4135
4136	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval_sp
4137	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'time_avg', &
4138	TRIM( time_average_text ) )
4139	ELSE
4140	nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
4141	TRIM( run_description_header ) )
4142	ENDIF
4143	CALL handle_netcdf_error( 'netcdf', 309 )
4144	nc_stat = NF90_ENDDEF( id_set_sp )
4145	CALL handle_netcdf_error( 'netcdf', 442 )
4146	message_string = 'NetCDF file for spectra ' // &
4147	'from previous run found.' // &
4148	'&This file will be extended.'
4149	CALL message( 'define_netcdf_header', 'PA0264', 0, 0, 0, 6, 0 )
4150
4151
4152	CASE ( 'pt_new' )
4153
4154	!
4155	!-- Define some global attributes of the dataset
4156	nc_stat = NF90_PUT_ATT( id_set_prt, NF90_GLOBAL, 'title', &
4157	TRIM( run_description_header ) )
4158	CALL handle_netcdf_error( 'netcdf', 310 )
4159
4160	!
4161	!-- Define time coordinate for particles (unlimited dimension)
4162	nc_stat = NF90_DEF_DIM( id_set_prt, 'time', NF90_UNLIMITED, &
4163	id_dim_time_prt )
4164	CALL handle_netcdf_error( 'netcdf', 311 )
4165
4166	nc_stat = NF90_DEF_VAR( id_set_prt, 'time', NF90_DOUBLE, &
4167	id_dim_time_prt, id_var_time_prt )
4168	CALL handle_netcdf_error( 'netcdf', 312 )
4169
4170	nc_stat = NF90_PUT_ATT( id_set_prt, id_var_time_prt, 'units', &
4171	'seconds' )
4172	CALL handle_netcdf_error( 'netcdf', 313 )
4173
4174	!
4175	!-- Define particle coordinate (maximum particle number)
4176	nc_stat = NF90_DEF_DIM( id_set_prt, 'particle_number', &
4177	maximum_number_of_particles, id_dim_prtnum )
4178	CALL handle_netcdf_error( 'netcdf', 314 )
4179
4180	nc_stat = NF90_DEF_VAR( id_set_prt, 'particle_number', NF90_DOUBLE, &
4181	id_dim_prtnum, id_var_prtnum )
4182	CALL handle_netcdf_error( 'netcdf', 315 )
4183
4184	nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prtnum, 'units', &
4185	'particle number' )
4186	CALL handle_netcdf_error( 'netcdf', 316 )
4187
4188	!
4189	!-- Define variable which contains the real number of particles in use
4190	nc_stat = NF90_DEF_VAR( id_set_prt, 'real_num_of_prt', NF90_INT, &
4191	id_dim_time_prt, id_var_rnop_prt )
4192	CALL handle_netcdf_error( 'netcdf', 317 )
4193
4194	nc_stat = NF90_PUT_ATT( id_set_prt, id_var_rnop_prt, 'units', &
4195	'particle number' )
4196	CALL handle_netcdf_error( 'netcdf', 318 )
4197
4198	!
4199	!-- Define the variables
4200	DO i = 1, 17
4201
4202	nc_stat = NF90_DEF_VAR( id_set_prt, prt_var_names(i), &
4203	nc_precision(8), &
4204	(/ id_dim_prtnum, id_dim_time_prt /), &
4205	id_var_prt(i) )
4206	CALL handle_netcdf_error( 'netcdf', 319 )
4207
4208	nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prt(i), &
4209	'long_name', TRIM( prt_var_names(i) ) )
4210	CALL handle_netcdf_error( 'netcdf', 320 )
4211
4212	nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prt(i), &
4213	'units', TRIM( prt_var_units(i) ) )
4214	CALL handle_netcdf_error( 'netcdf', 321 )
4215
4216	ENDDO
4217
4218	!
4219	!-- Leave NetCDF define mode
4220	nc_stat = NF90_ENDDEF( id_set_prt )
4221	CALL handle_netcdf_error( 'netcdf', 322 )
4222
4223
4224	CASE ( 'pt_ext' )
4225
4226	!
4227	!-- Get the id of the time coordinate (unlimited coordinate) and its
4228	!-- last index on the file. The next time level is prt..count+1.
4229	!-- The current time must be larger than the last output time
4230	!-- on the file.
4231	nc_stat = NF90_INQ_VARID( id_set_prt, 'time', id_var_time_prt )
4232	CALL handle_netcdf_error( 'netcdf', 323 )
4233
4234	nc_stat = NF90_INQUIRE_VARIABLE( id_set_prt, id_var_time_prt, &
4235	dimids = id_dim_time_old )
4236	CALL handle_netcdf_error( 'netcdf', 324 )
4237	id_dim_time_prt = id_dim_time_old(1)
4238
4239	nc_stat = NF90_INQUIRE_DIMENSION( id_set_prt, id_dim_time_prt, &
4240	len = prt_time_count )
4241	CALL handle_netcdf_error( 'netcdf', 325 )
4242
4243	nc_stat = NF90_GET_VAR( id_set_prt, id_var_time_prt, &
4244	last_time_coordinate, &
4245	start = (/ prt_time_count /), &
4246	count = (/ 1 /) )
4247	CALL handle_netcdf_error( 'netcdf', 326 )
4248
4249	IF ( last_time_coordinate(1) >= simulated_time ) THEN
4250	message_string = 'NetCDF file for particles ' // &
4251	'from previous run found,' // &
4252	'&but this file cannot be extended becaus' // &
4253	'e the current output time' // &
4254	'&is less or equal than the last output t' // &
4255	'ime on this file.' // &
4256	'&New file is created instead.'
4257	CALL message( 'define_netcdf_header', 'PA0265', 0, 1, 0, 6, 0 )
4258	prt_time_count = 0
4259	extend = .FALSE.
4260	RETURN
4261	ENDIF
4262
4263	!
4264	!-- Dataset seems to be extendable.
4265	!-- Now get the variable ids.
4266	nc_stat = NF90_INQ_VARID( id_set_prt, 'real_num_of_prt', &
4267	id_var_rnop_prt )
4268	CALL handle_netcdf_error( 'netcdf', 327 )
4269
4270	DO i = 1, 17
4271
4272	nc_stat = NF90_INQ_VARID( id_set_prt, prt_var_names(i), &
4273	id_var_prt(i) )
4274	CALL handle_netcdf_error( 'netcdf', 328 )
4275
4276	ENDDO
4277
4278	message_string = 'NetCDF file for particles ' // &
4279	'from previous run found.' // &
4280	'&This file will be extended.'
4281	CALL message( 'define_netcdf_header', 'PA0266', 0, 0, 0, 6, 0 )
4282
4283
4284
4285	CASE ( 'ps_new' )
4286
4287	!
4288	!-- Define some global attributes of the dataset
4289	nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
4290	TRIM( run_description_header ) )
4291	CALL handle_netcdf_error( 'netcdf', 396 )
4292
4293	!
4294	!-- Define time coordinate for particle time series (unlimited dimension)
4295	nc_stat = NF90_DEF_DIM( id_set_pts, 'time', NF90_UNLIMITED, &
4296	id_dim_time_pts )
4297	CALL handle_netcdf_error( 'netcdf', 397 )
4298
4299	nc_stat = NF90_DEF_VAR( id_set_pts, 'time', NF90_DOUBLE, &
4300	id_dim_time_pts, id_var_time_pts )
4301	CALL handle_netcdf_error( 'netcdf', 398 )
4302
4303	nc_stat = NF90_PUT_ATT( id_set_pts, id_var_time_pts, 'units', &
4304	'seconds')
4305	CALL handle_netcdf_error( 'netcdf', 399 )
4306
4307	!
4308	!-- Define the variables. If more than one particle group is defined,
4309	!-- define seperate variables for each group
4310	var_list = ';'
4311	DO i = 1, dopts_num
4312
4313	!
4314	!-- First, remove those characters not allowed by NetCDF
4315	netcdf_var_name = dopts_label(i)
4316	CALL clean_netcdf_varname( netcdf_var_name )
4317
4318	DO j = 0, number_of_particle_groups
4319
4320	IF ( j == 0 ) THEN
4321	suffix1 = ''
4322	ELSE
4323	WRITE ( suffix1, '(''_'',I2.2)' ) j
4324	ENDIF
4325
4326	nc_stat = NF90_DEF_VAR( id_set_pts, &
4327	TRIM( netcdf_var_name ) // suffix1, &
4328	nc_precision(6), &
4329	(/ id_dim_time_pts /), &
4330	id_var_dopts(i,j) )
4331	CALL handle_netcdf_error( 'netcdf', 400 )
4332
4333	IF ( j == 0 ) THEN
4334	nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
4335	'long_name', &
4336	TRIM( dopts_label(i) ) )
4337	ELSE
4338	nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
4339	'long_name', &
4340	TRIM( dopts_label(i) ) // ' PG ' // &
4341	suffix1(2:3) )
4342	ENDIF
4343	CALL handle_netcdf_error( 'netcdf', 401 )
4344
4345	nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
4346	'units', TRIM( dopts_unit(i) ) )
4347	CALL handle_netcdf_error( 'netcdf', 402 )
4348
4349	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
4350	suffix1 // ';'
4351
4352	IF ( number_of_particle_groups == 1 ) EXIT
4353
4354	ENDDO
4355
4356	ENDDO
4357
4358	!
4359	!-- Write the list of variables as global attribute (this is used by
4360	!-- restart runs)
4361	nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
4362	var_list )
4363	CALL handle_netcdf_error( 'netcdf', 403 )
4364
4365
4366	!
4367	!-- Leave NetCDF define mode
4368	nc_stat = NF90_ENDDEF( id_set_pts )
4369	CALL handle_netcdf_error( 'netcdf', 404 )
4370
4371
4372	CASE ( 'ps_ext' )
4373
4374	!
4375	!-- Get the list of variables and compare with the actual run.
4376	!-- First var_list_old has to be reset, since GET_ATT does not assign
4377	!-- trailing blanks.
4378	var_list_old = ' '
4379	nc_stat = NF90_GET_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
4380	var_list_old )
4381	CALL handle_netcdf_error( 'netcdf', 405 )
4382
4383	var_list = ';'
4384	i = 1
4385	DO i = 1, dopts_num
4386
4387	netcdf_var_name = dopts_label(i)
4388	CALL clean_netcdf_varname( netcdf_var_name )
4389
4390	DO j = 0, number_of_particle_groups
4391
4392	IF ( j == 0 ) THEN
4393	suffix1 = ''
4394	ELSE
4395	WRITE ( suffix1, '(''_'',I2.2)' ) j
4396	ENDIF
4397
4398	var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
4399	suffix1 // ';'
4400
4401	IF ( number_of_particle_groups == 1 ) EXIT
4402
4403	ENDDO
4404
4405	ENDDO
4406
4407	IF ( TRIM( var_list ) /= TRIM( var_list_old ) ) THEN
4408	message_string = 'NetCDF file for particle time series ' // &
4409	'from previous run found,' // &
4410	'& but this file cannot be extended due to' // &
4411	' variable mismatch.' // &
4412	'&New file is created instead.'
4413	CALL message( 'define_netcdf_header', 'PA0267', 0, 1, 0, 6, 0 )
4414	extend = .FALSE.
4415	RETURN
4416	ENDIF
4417
4418	!
4419	!-- Get the id of the time coordinate (unlimited coordinate) and its
4420	!-- last index on the file. The next time level is dots..count+1.
4421	!-- The current time must be larger than the last output time
4422	!-- on the file.
4423	nc_stat = NF90_INQ_VARID( id_set_pts, 'time', id_var_time_pts )
4424	CALL handle_netcdf_error( 'netcdf', 406 )
4425
4426	nc_stat = NF90_INQUIRE_VARIABLE( id_set_pts, id_var_time_pts, &
4427	dimids = id_dim_time_old )
4428	CALL handle_netcdf_error( 'netcdf', 407 )
4429	id_dim_time_pts = id_dim_time_old(1)
4430
4431	nc_stat = NF90_INQUIRE_DIMENSION( id_set_pts, id_dim_time_pts, &
4432	len = dopts_time_count )
4433	CALL handle_netcdf_error( 'netcdf', 408 )
4434
4435	nc_stat = NF90_GET_VAR( id_set_pts, id_var_time_pts, &
4436	last_time_coordinate, &
4437	start = (/ dopts_time_count /), &
4438	count = (/ 1 /) )
4439	CALL handle_netcdf_error( 'netcdf', 409 )
4440
4441	IF ( last_time_coordinate(1) >= simulated_time ) THEN
4442	message_string = 'NetCDF file for particle time series ' // &
4443	'from previous run found,' // &
4444	'&but this file cannot be extended becaus' // &
4445	'e the current output time' // &
4446	'&is less or equal than the last output t' // &
4447	'ime on this file.' // &
4448	'&New file is created instead.'
4449	CALL message( 'define_netcdf_header', 'PA0268', 0, 1, 0, 6, 0 )
4450	dopts_time_count = 0
4451	extend = .FALSE.
4452	RETURN
4453	ENDIF
4454
4455	!
4456	!-- Dataset seems to be extendable.
4457	!-- Now get the variable ids
4458	i = 1
4459	DO i = 1, dopts_num
4460
4461	netcdf_var_name_base = dopts_label(i)
4462	CALL clean_netcdf_varname( netcdf_var_name_base )
4463
4464	DO j = 0, number_of_particle_groups
4465
4466	IF ( j == 0 ) THEN
4467	suffix1 = ''
4468	ELSE
4469	WRITE ( suffix1, '(''_'',I2.2)' ) j
4470	ENDIF
4471
4472	netcdf_var_name = TRIM( netcdf_var_name_base ) // suffix1
4473
4474	nc_stat = NF90_INQ_VARID( id_set_pts, netcdf_var_name, &
4475	id_var_dopts(i,j) )
4476	CALL handle_netcdf_error( 'netcdf', 410 )
4477
4478	IF ( number_of_particle_groups == 1 ) EXIT
4479
4480	ENDDO
4481
4482	ENDDO
4483
4484	!
4485	!-- Update the title attribute on file
4486	!-- In order to avoid 'data mode' errors if updated attributes are larger
4487	!-- than their original size, NF90_PUT_ATT is called in 'define mode'
4488	!-- enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
4489	!-- performance loss due to data copying; an alternative strategy would be
4490	!-- to ensure equal attribute size in a job chain. Maybe revise later.
4491	nc_stat = NF90_REDEF( id_set_pts )
4492	CALL handle_netcdf_error( 'netcdf', 443 )
4493	nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
4494	TRIM( run_description_header ) )
4495	CALL handle_netcdf_error( 'netcdf', 411 )
4496	nc_stat = NF90_ENDDEF( id_set_pts )
4497	CALL handle_netcdf_error( 'netcdf', 444 )
4498	message_string = 'NetCDF file for particle time series ' // &
4499	'from previous run found.' // &
4500	'&This file will be extended.'
4501	CALL message( 'define_netcdf_header', 'PA0269', 0, 0, 0, 6, 0 )
4502
4503
4504	CASE DEFAULT
4505
4506	message_string = 'mode "' // TRIM( mode) // '" not supported'
4507	CALL message( 'define_netcdf_header', 'PA0270', 0, 0, 0, 6, 0 )
4508
4509	END SELECT
4510
4511	#endif
4512	END SUBROUTINE define_netcdf_header
4513
4514
4515	SUBROUTINE handle_netcdf_error( routine_name, errno )
4516	#if defined( __netcdf )
4517
4518	!------------------------------------------------------------------------------!
4519	!
4520	! Description:
4521	! ------------
4522	! Prints out a text message corresponding to the current status.
4523	!------------------------------------------------------------------------------!
4524
4525	USE control_parameters
4526	USE netcdf
4527	USE netcdf_control
4528	USE pegrid
4529
4530	IMPLICIT NONE
4531
4532	CHARACTER(LEN=6) :: message_identifier
4533	CHARACTER(LEN=*) :: routine_name
4534
4535	INTEGER :: errno
4536
4537	IF ( nc_stat /= NF90_NOERR ) THEN
4538
4539	WRITE( message_identifier, '(''NC'',I4.4)' ) errno
4540	message_string = TRIM( NF90_STRERROR( nc_stat ) )
4541
4542	CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
4543
4544	ENDIF
4545
4546	#endif
4547	END SUBROUTINE handle_netcdf_error
4548
4549
4550
4551	SUBROUTINE clean_netcdf_varname( string )
4552	#if defined( __netcdf )
4553
4554	!------------------------------------------------------------------------------!
4555	!
4556	! Description:
4557	! ------------
4558	! Replace those characters in string which are not allowed by NetCDF.
4559	!------------------------------------------------------------------------------!
4560
4561	USE netcdf_control
4562
4563	IMPLICIT NONE
4564
4565	CHARACTER (LEN=10), INTENT(INOUT) :: string
4566
4567	INTEGER :: i, ic
4568
4569	DO i = 1, replace_num
4570	DO
4571	ic = INDEX( string, replace_char(i) )
4572	IF ( ic /= 0 ) THEN
4573	string(ic:ic) = replace_by(i)
4574	ELSE
4575	EXIT
4576	ENDIF
4577	ENDDO
4578	ENDDO
4579
4580	#endif
4581	END SUBROUTINE clean_netcdf_varname

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |