Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |