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

Last change on this file since 2040 was 2040, checked in by gronemeier, 5 years ago

last commit documented

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