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

Last change on this file since 2038 was 2038, checked in by knoop, 5 years ago

last commit documented

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