source: palm/trunk/SOURCE/netcdf_interface.f90 @ 1786

Last change on this file since 1786 was 1786, checked in by raasch, 8 years ago

pmc-change in server-client get-put, spectra-directives removed, spectra-package modularized

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