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

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

last commit documented

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