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

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

last commit documented

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