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

Last change on this file since 3582 was 3582, checked in by suehring, 6 years ago

Merge branch salsa with trunk

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