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

Last change on this file since 4178 was 4127, checked in by suehring, 2 years ago

Merge with branch resler: biomet- output of bio_mrt added; plant_canopy - separate vertical dimension for 3D output (to save disk space); radiation - remove unused plant canopy variables; urban-surface model - do not add anthropogenic heat during wall spin-up

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