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

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

added _mod string to several filenames to meet the naming convection for modules

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