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

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

changes related to steering and formating of urban surface model

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