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

Last change on this file since 4029 was 4029, checked in by raasch, 2 years ago

bugfix: decycling of chemistry species after nesting data transfer, exchange of ghost points and boundary conditions separated for chemical species and SALSA module, nest_chemistry option removed, netcdf variable NF90_NOFILL is used as argument instead of 1 in calls to NF90_DEF_VAR_FILL

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