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

Last change on this file since 3187 was 3187, checked in by sward, 3 years ago

Reworked agent pathfinding and modified output

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