Home

Context Navigation

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

Last change on this file since 1745 was 1745, checked in by gronemeier, 8 years ago
Bugfix:calculation of time levels for parallel NetCDF output
Property svn:keywords set to `Id`
File size: 219.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |