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

Last change on this file since 2012 was 2012, checked in by kanani, 5 years ago

last commit documented

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