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

Last change on this file since 3459 was 3459, checked in by gronemeier, 3 years ago

add longitude and latitude to output files; update run_control files for tests

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