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

Last change on this file since 3942 was 3942, checked in by kanani, 2 years ago

Fix too short driver attribute lengths, and individualize error messages (netcdf_data_input_mod, netcdf_interface_mod)

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