Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |