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

Last change on this file since 3294 was 3294, checked in by raasch, 3 years ago

modularization of the ocean code

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