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

Last change on this file since 2746 was 2746, checked in by suehring, 3 years ago

Read information from statitic driver for resolved vegetation independently from land- or urban-surface model

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