Home

Context Navigation

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

Last change on this file since 1295 was 1207, checked in by witha, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 190.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |