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

Last change on this file since 3235 was 3235, checked in by sward, 3 years ago

Bugfix in MAS output, added related messages, reworked MAS cpu logging

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