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

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

diagnostic output: Modularize diagnostic output, rename subroutines; formatting adjustments; allocate arrays only when required; add output of uu, vv, ww to enable variance calculation via temporal EC method; radiation: bugfix in masked data output; flow_statistics: Correct conversion to kinematic vertical scalar fluxes in case of pw-scheme and statistic regions

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