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

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

Added agent ID output

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