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

Last change on this file since 3995 was 3995, checked in by suehring, 2 years ago

remove unused variable and fix non-standard string operation PGI compiler

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