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

Last change on this file since 2512 was 2512, checked in by raasch, 7 years ago

upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny; no output if redundant ghost layer data to NetCDF files

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