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

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

further modularization of land surface model (2D/3D output and restart data). Bugfix for restart runs without land surface model

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