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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

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