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

Last change on this file since 2696 was 2696, checked in by kanani, 4 years ago

Merge of branch palm4u into trunk

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