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

Last change on this file since 2210 was 2210, checked in by kanani, 4 years ago

last commit documented

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