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

Last change on this file since 2039 was 2039, checked in by gronemeier, 5 years ago

Increased the number of possible statistic regions to 99

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