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

Last change on this file since 3768 was 3766, checked in by raasch, 2 years ago

unused_variables removed, bugfix in im_define_netcdf_grid argument list, trim added to avoid truncation compiler warnings, save attribute added to local targets to avoid outlive pointer target warning, first argument removed from module_interface_rrd_*, file module_interface reformatted with respect to coding standards, bugfix in surface_data_output_rrd_local (variable k removed)

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