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

Last change on this file since 1783 was 1783, checked in by raasch, 5 years ago

NetCDF routines modularized; new parameter netcdf_deflate; further changes in the pmc

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