Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |