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

Last change on this file since 2037 was 2037, checked in by knoop, 5 years ago

Anelastic approximation implemented

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