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

Last change on this file since 2716 was 2716, checked in by kanani, 3 years ago

Correction of "Former revisions" section

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