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

Last change on this file since 1972 was 1972, checked in by maronga, 5 years ago

further modularization of land surface model (2D/3D output and restart data)

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