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

Last change on this file since 4069 was 4069, checked in by Giersch, 2 years ago

Bugfix for masked output, compiler warning removed, test case for wind turbine model revised

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