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

Last change on this file since 2292 was 2292, checked in by schwenkel, 4 years ago

implementation of new bulk microphysics scheme

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