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

Last change on this file since 2769 was 2769, checked in by raasch, 3 years ago

bugfix for calculating number of required output time levels in case of output at the beginning of a restart run

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