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

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

Reading of 3D topography using NetCDF data type NC_BYTE; bugfixes in reading 3D topography from file

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