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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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