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

Last change on this file since 4196 was 4196, checked in by gronemeier, 22 months ago

Consider rotation of model domain for claculation of Coriolis force:

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