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

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

bugfix: initialize return values to ensure they are set before returning (routine define_geo_coordinates); change order of dimensions for some variables

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