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

Last change on this file since 3045 was 3045, checked in by Giersch, 3 years ago

Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled

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