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

Last change on this file since 3467 was 3467, checked in by suehring, 3 years ago

Branch salsa @3446 re-integrated into trunk

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