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

Last change on this file since 3421 was 3421, checked in by gronemeier, 3 years ago

new surface-data output; renamed output variables (pt to theta, rho_air to rho, rho_ocean to rho_sea_water)

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