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

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

flight module added

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