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

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

Coupling of indoor model to atmosphere; output of indoor temperatures and waste heat; enable restarts with indoor model; bugfix plant transpiration; bugfix - missing calculation of 10cm temperature

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