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

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

last commit documented

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