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

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

small bugfix, formatting and new PCM output

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