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

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

last commit documented

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