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

Last change on this file since 3198 was 3198, checked in by sward, 6 years ago

Added MAS end time, used time_since_reference_point, corrected tolerance_dp in nav_mesh

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