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

Last change on this file since 1987 was 1981, checked in by suehring, 8 years ago

last commit documented

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