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

Last change on this file since 3727 was 3727, checked in by gronemeier, 2 years ago

New: surface output available in NetCDF format (Makefile, netcdf_interface_mod, surface_data_output_mod)

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