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

Last change on this file since 2964 was 2964, checked in by Giersch, 3 years ago

Bugfix in the calculation of fixed number of output time levels for parallel netcdf output, error message related to reading sky view factors revised

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