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

Last change on this file since 2265 was 2265, checked in by schwenkel, 4 years ago

unused variables removed

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