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

Last change on this file since 2270 was 2270, checked in by maronga, 4 years ago

major revisions in land surface model

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