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

Last change on this file since 4046 was 4046, checked in by knoop, 2 years ago

removal of special treatment for usm_define_netcdf_grid call in netcdf_interface_mod

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