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

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

from branch resler@3462: add MRT shaping function (radiation_model_mod), use basic constants (biometeorology_mod), adjust precision to wp (biometeorology_pt_mod)

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