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

Last change on this file since 2032 was 2032, checked in by knoop, 5 years ago

last commit documented

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