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

Last change on this file since 2007 was 2007, checked in by kanani, 5 years ago

changes in the course of urban surface model implementation

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