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

Last change on this file since 3525 was 3525, checked in by kanani, 3 years ago

Changes related to clean-up of biometeorology (by dom_dwd_user)

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