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

Last change on this file since 3337 was 3337, checked in by kanani, 3 years ago

reintegrate branch resler to trunk

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