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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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