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

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

last commit documented

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