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

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

Fix for biomet output (ticket:757), merge of uv_exposure into biometeorology_mod

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