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

Last change on this file since 2932 was 2932, checked in by maronga, 3 years ago

renamed all Fortran NAMELISTS

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