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

Last change on this file since 3665 was 3665, checked in by raasch, 2 years ago

dummy statements added to avoid compiler warnings about unused variables, unused variables removed, ssh-call for submitting batch jobs on remote systems modified again to avoid output of login messages on specific systems

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