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

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

last commit documented

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