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

Last change on this file since 2008 was 2008, checked in by kanani, 5 years ago

last commit documented

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