source: palm/trunk/SOURCE/netcdf_interface_mod.f90

Last change on this file was 4856, checked in by raasch, 5 months ago

array sizes for output profiles increased from 400 to 500

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