Home

Context Navigation

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

Last change on this file since 1691 was 1691, checked in by maronga, 10 years ago
various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings
Property svn:keywords set to `Id`
File size: 212.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |