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

Last change on this file since 2109 was 2109, checked in by raasch, 4 years ago

bugfix: length of character string netcdf_var_name extended to avoid problems with restart runs due to truncation

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