source: palm/trunk/SOURCE/netcdf_interface.f90 @ 1818

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

last commit documented / copyright update

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