Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |