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

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

various changes to avoid compiler warnings (mainly removal of unused variables)

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