Home

Context Navigation

netcdf.f90 @ 493

Last change on this file since 493 was 493, checked in by raasch, 14 years ago

New:
---
Output in NetCDF4-format. New d3par-parameter netcdf_data_format.

(check_open, check_parameters, close_file, data_output_2d, data_output_3d, header, modules, netcdf, parin)

Modules to be loaded for compilation (mbuild) or job execution (mrun)
can be given in the configuration file using variable modules. Example:

%modules ifort/11.0.069:netcdf lcsgih parallel

This method replaces the (undocumented) mpilib-variable.

WARNING: All fixed settings of modules in the scripts mbuild, mrun, and subjob
have been removed! Please set the modules variable appropriately in your
configuration file. (mbuild, mrun, subjob)

Changed:

Parameters netcdf_64bit and netcdf_64bit_3d have been removed. Use
netcdf_data_format = 2 for choosing the classic 64bit-offset format (this is
the default). The offset-format can not be set independently for the
3d-output-data any more.

Parameters netcdf_format_mask, netcdf_format_mask_av, and variables
nc_format_mask, format_parallel_io removed. They are replaced by the new
parameter netcdf_data_format. (check_open, close_file,
data_output_mask, header, init_masks, modules, parin)

Errors:

bugfix in trunk/UTIL/Makefile: forgot to compile for interpret_config

Bugfix: timeseries data have to be collected by PE0 (user_statistics)

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

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

Context Navigation

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

Download in other formats: