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

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

last commit documented

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