Home

Context Navigation

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

Last change on this file since 1070 was 1054, checked in by hoffmann, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 190.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |