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

Last change on this file since 2110 was 2110, checked in by raasch, 4 years ago

last commit documented

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