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

Last change on this file since 3159 was 3159, checked in by sward, 7 years ago

Added multi agent system

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