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

Last change on this file since 2239 was 2239, checked in by suehring, 4 years ago

Bugfix xy-output of land-surface variables

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