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

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

Separate balance equations for humidity and passive_scalar

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