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

Last change on this file since 3745 was 3745, checked in by suehring, 2 years ago

document last commit

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