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

Last change on this file since 3004 was 3004, checked in by Giersch, 3 years ago

precipitation_rate removed, further allocation checks for data output of averaged quantities implemented, double CALL of flow_statistics at the beginning of time_integration removed, further minor bugfixes, comments added

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