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

Last change on this file since 2766 was 2766, checked in by kanani, 3 years ago

Removal of chem directive, plus minor changes

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