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

Last change on this file since 2817 was 2817, checked in by knoop, 3 years ago

Preliminary gust module interface implemented

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