source: palm/trunk/SOURCE/netcdf_interface_mod.f90 @ 1980

Last change on this file since 1980 was 1980, checked in by suehring, 5 years ago

Bugfix, in order to steer user-defined output, setting flag found explicitly

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