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

Last change on this file since 1961 was 1961, checked in by suehring, 5 years ago

last commit documented

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