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

Last change on this file since 1991 was 1991, checked in by gronemeier, 5 years ago

last commit documented

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