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

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

Moved user_define_netdf_grid into user_module.f90
Added module interface for the definition of additional timeseries

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