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

Last change on this file since 4598 was 4502, checked in by schwenkel, 12 months ago

Implementation of ice microphysics

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