Home

Context Navigation

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

Last change on this file since 1053 was 1053, checked in by hoffmann, 11 years ago
two-moment cloud physics implemented
Property svn:keywords set to `Id`
File size: 190.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |