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

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

Revision history corrected

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