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

Last change on this file since 1850 was 1850, checked in by maronga, 5 years ago

added _mod string to several filenames to meet the naming convection for modules

  • Property svn:keywords set to Id
File size: 218.9 KB
RevLine 
[1850]1!> @file netcdf_interface_mod.f90
[1036]2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
[1818]16! Copyright 1997-2016 Leibniz Universitaet Hannover
[1036]17!--------------------------------------------------------------------------------!
18!
[1]19! Current revisions:
20! ------------------
[1850]21! Module renamed
[1834]22!
23!
[1784]24! Former revisions:
25! -----------------
26! $Id: netcdf_interface_mod.f90 1850 2016-04-08 13:29:27Z maronga $
27!
[1834]28! 1833 2016-04-07 14:23:03Z raasch
29! spectrum renamed spectra_mod
30!
[1787]31! 1786 2016-03-08 05:49:27Z raasch
32! Bugfix: id_var_time_sp made public
33!
[1784]34! 1783 2016-03-06 18:36:17Z raasch
[1783]35! netcdf interface has been modularized, former file netcdf renamed to
36! netcdf_interface, creation of netcdf-dimensions and -variables moved to
37! specific new subroutines create_netcdf_dim and create_netcdf_var,
38! compression (deflation) of variables implemented,
39! ibmy special cpp directive removed
[1321]40!
[1746]41! 1745 2016-02-05 13:06:51Z gronemeier
42! Bugfix: recalculating ntdim_3d, ntdim_2d_xy/xz/yz when checking the
43!         extensibility of an existing file (only when using parallel NetCDF).
44!
[1692]45! 1691 2015-10-26 16:17:44Z maronga
46! Added output of radiative heating rates for RRTMG. Corrected output of
47! radiative fluxes
48!
[1683]49! 1682 2015-10-07 23:56:08Z knoop
50! Code annotations made doxygen readable
51!
[1597]52! 1596 2015-05-21 09:34:28Z gronemeier
53! Bugfix in masked data output. Read 'zu_3d' when trying to extend masked data
54!
[1552]55! 1551 2015-03-03 14:18:16Z maronga
56! Added support for land surface model and radiation model output. In the course
57! of this action a new vertical grid zs (soil) was introduced.
58!
[1354]59! 1353 2014-04-08 15:21:23Z heinze
60! REAL constants provided with KIND-attribute
61!
[1323]62! 1322 2014-03-20 16:38:49Z raasch
63! Forgotten ONLY-attribute added to USE-statements
64!
[1321]65! 1320 2014-03-20 08:40:49Z raasch
[1320]66! ONLY-attribute added to USE-statements,
67! kind-parameters added to all INTEGER and REAL declaration statements,
68! kinds are defined in new module kinds,
69! revision history before 2012 removed,
70! comment fields (!:) to be used for variable explanations added to
71! all variable declaration statements
[1309]72!
73! 1308 2014-03-13 14:58:42Z fricke
[1308]74! +ntime_count, oldmode
75! Adjust NF90_CREATE and NF90_OPEN statement for parallel output
76! To increase the performance for parallel output, the following is done:
77! - Limit time dimension
78! - Values of axis data are only written by PE0
79! - No fill is set for all variables
80! Check the number of output time levels for restart jobs
[520]81!
[1207]82! 1206 2013-07-18 12:49:16Z witha
83! Bugfix: typo in preprocessor directive in subroutine open_write_netcdf_file
84!
[1093]85! 1092 2013-02-02 11:24:22Z raasch
86! unused variables removed
87!
[1054]88! 1053 2012-11-13 17:11:03Z hoffmann
89! +qr, nr, prr
90!
[1037]91! 1036 2012-10-22 13:43:42Z raasch
92! code put under GPL (PALM 3.9)
93!
[1035]94! 1031 2012-10-19 14:35:30Z raasch
95! netCDF4 without parallel file support implemented, new routines
96! create_netcdf_file and open_write_netcdf_file at end
97!
[993]98! 992 2012-09-05 15:08:26Z hoffmann
99! Removal of the informative messages PA0352 and PA0353.
[984]100
101! 983 2012-08-21 14:17:57Z hoffmann
102! Bugfix in cross_profiles.
[520]103!
[965]104! 964 2012-07-26 09:14:24Z raasch
105! rev 951 and 959 reformatted
106!
[960]107! 959 2012-07-24 13:13:41Z hoffmann
[964]108! Bugfix in cross_profiles. It is not allowed to arrange more than 100
[960]109! profiles with cross_profiles.
110!
[952]111! 951 2012-07-19 14:22:52Z hoffmann
[1031]112! cross_profiles, profile_rows, profile_columns are written to netCDF header
[952]113!
[1]114! Revision 1.1  2005/05/18 15:37:16  raasch
115! Initial revision
116!
117!
118! Description:
119! ------------
[1682]120!> In case of extend = .FALSE.:
121!> Define all necessary dimensions, axes and variables for the different
122!> netCDF datasets. This subroutine is called from check_open after a new
123!> dataset is created. It leaves the open netCDF files ready to write.
124!>
125!> In case of extend = .TRUE.:
126!> Find out if dimensions and variables of an existing file match the values
127!> of the actual run. If so, get all necessary information (ids, etc.) from
128!> this file.
129!>
130!> Parameter av can assume values 0 (non-averaged data) and 1 (time averaged
131!> data)
[1745]132!>
133!> @todo calculation of output time levels for parallel NetCDF still does not
134!>       cover every exception (change of dt_do, end_time in restart)
[1]135!------------------------------------------------------------------------------!
[1783]136 MODULE netcdf_interface
137
138    USE control_parameters, ONLY: max_masks
139    USE kinds
140#if defined( __netcdf )
141    USE NETCDF
[1682]142#endif
[1783]143
144    PRIVATE
145
146    INTEGER(iwp), PARAMETER ::  dopr_norm_num = 7, dopts_num = 29, dots_max = 100
147
148    CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_names =          &
149         (/ 'wpt0  ', 'ws2   ', 'tsw2  ', 'ws3   ', 'ws2tsw', 'wstsw2',        &
150            'z_i   ' /)
151
152    CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_longnames =      &
153         (/ 'wpt0  ', 'w*2   ', 't*w2  ', 'w*3   ', 'w*2t*w', 'w*t*w2',        &
154            'z_i   ' /)
155
156    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_label =                   &
157          (/ 'tnpt   ', 'x_     ', 'y_     ', 'z_     ', 'z_abs  ', 'u      ', &
158             'v      ', 'w      ', 'u"     ', 'v"     ', 'w"     ', 'npt_up ', &
159             'w_up   ', 'w_down ', 'radius ', 'r_min  ', 'r_max  ', 'npt_max', &
160             'npt_min', 'x*2    ', 'y*2    ', 'z*2    ', 'u*2    ', 'v*2    ', &
161             'w*2    ', 'u"2    ', 'v"2    ', 'w"2    ', 'npt*2  ' /)
162
163    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_unit =                    &
164          (/ 'number ', 'm      ', 'm      ', 'm      ', 'm      ', 'm/s    ', &
165             'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'number ', &
166             'm/s    ', 'm/s    ', 'm      ', 'm      ', 'm      ', 'number ', &
167             'number ', 'm2     ', 'm2     ', 'm2     ', 'm2/s2  ', 'm2/s2  ', &
168             'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'number2' /)
169
170    INTEGER(iwp) ::  dots_num  = 31  !< number of timeseries defined by default
171    INTEGER(iwp) ::  dots_soil = 24  !< starting index for soil-timeseries
172    INTEGER(iwp) ::  dots_rad  = 32  !< starting index for radiation-timeseries
173
174    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label =                    &
175          (/ 'E            ', 'E*           ', 'dt           ',                &
176             'u*           ', 'th*          ', 'umax         ',                &
177             'vmax         ', 'wmax         ', 'div_new      ',                &
178             'div_old      ', 'z_i_wpt      ', 'z_i_pt       ',                &
179             'w*           ', 'w"pt"0       ', 'w"pt"        ',                &
180             'wpt          ', 'pt(0)        ', 'pt(z_mo)     ',                &
181             'w"u"0        ', 'w"v"0        ', 'w"q"0        ',                &
182             'ol           ', 'q*           ', 'ghf_eb       ',                &
183             'shf_eb       ', 'qsws_eb      ', 'qsws_liq_eb  ',                &
184             'qsws_soil_eb ', 'qsws_veg_eb  ', 'r_a          ',                &
185             'r_s          ', 'rad_net      ', 'rad_lw_in    ',                &
186             'rad_lw_out   ', 'rad_sw_in    ', 'rad_sw_out   ',                &
187             'rrtm_aldif   ', 'rrtm_aldir   ', 'rrtm_asdif   ',                &
188             'rrtm_asdir   ',                                                  &
189             ( 'unknown      ', i9 = 1, dots_max-40 ) /)
190
191    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit =                     &
192          (/ 'm2/s2        ', 'm2/s2        ', 's            ',                &
193             'm/s          ', 'K            ', 'm/s          ',                &
194             'm/s          ', 'm/s          ', 's-1          ',                &
195             's-1          ', 'm            ', 'm            ',                &
196             'm/s          ', 'K m/s        ', 'K m/s        ',                &
197             'K m/s        ', 'K            ', 'K            ',                &
198             'm2/s2        ', 'm2/s2        ', 'kg m/s       ',                &
199             'm            ', 'kg/kg        ', '             ',                &
200             '             ', '             ', '             ',                &
201             '             ', 'W/m2         ', 's/m          ',                &
202             '             ', 'W/m2         ', 'W/m2         ',                &
203             'W/m2         ', 'W/m2         ', 'W/m2         ',                &
204             '             ', '             ', '             ',                &
205             '             ',                                                  &
206             ( 'unknown      ', i9 = 1, dots_max-40 ) /)
207
208    CHARACTER (LEN=9), DIMENSION(300) ::  dopr_unit = 'unknown'
209
210    CHARACTER (LEN=7), DIMENSION(0:1,100) ::  do2d_unit, do3d_unit
211
212    CHARACTER (LEN=16), DIMENSION(25) ::  prt_var_names = &
213          (/ 'pt_age          ', 'pt_dvrp_size    ', 'pt_origin_x     ', &
214             'pt_origin_y     ', 'pt_origin_z     ', 'pt_radius       ', &
215             'pt_speed_x      ', 'pt_speed_y      ', 'pt_speed_z      ', &
216             'pt_weight_factor', 'pt_x            ', 'pt_y            ', &
217             'pt_z            ', 'pt_color        ', 'pt_group        ', &
218             'pt_tailpoints   ', 'pt_tail_id      ', 'pt_density_ratio', &
219             'pt_exp_arg      ', 'pt_exp_term     ', 'not_used        ', &
220             'not_used        ', 'not_used        ', 'not_used        ', &
221             'not_used        ' /)
222
223    CHARACTER (LEN=16), DIMENSION(25) ::  prt_var_units = &
224          (/ 'seconds         ', 'meters          ', 'meters          ', &
225             'meters          ', 'meters          ', 'meters          ', &
226             'm/s             ', 'm/s             ', 'm/s             ', &
227             'factor          ', 'meters          ', 'meters          ', &
228             'meters          ', 'none            ', 'none            ', &
229             'none            ', 'none            ', 'ratio           ', &
230             'none            ', 'none            ', 'not_used        ', &
231             'not_used        ', 'not_used        ', 'not_used        ', &
232             'not_used        ' /)
233
234    CHARACTER(LEN=20), DIMENSION(11) ::  netcdf_precision = ' '
235    CHARACTER(LEN=40) ::  netcdf_data_format_string
236
237    INTEGER(iwp) ::  id_dim_prtnum, id_dim_time_pr, id_dim_time_prt, &
238                     id_dim_time_pts, id_dim_time_sp, id_dim_time_ts, id_dim_x_sp, &
239                     id_dim_y_sp, id_dim_zu_sp, id_dim_zw_sp, id_set_pr, &
240                     id_set_prt, id_set_pts, id_set_sp, id_set_ts, id_var_prtnum, &
241                     id_var_rnop_prt, id_var_time_pr, id_var_time_prt, &
242                     id_var_time_pts, id_var_time_sp, id_var_time_ts, id_var_x_sp, &
243                     id_var_y_sp, id_var_zu_sp, id_var_zw_sp, nc_stat
244
245    INTEGER(iwp), DIMENSION(0:1) ::  id_dim_time_xy, id_dim_time_xz, &
246                    id_dim_time_yz, id_dim_time_3d, id_dim_x_xy, id_dim_xu_xy, &
247                    id_dim_x_xz, id_dim_xu_xz, id_dim_x_yz, id_dim_xu_yz, &
248                    id_dim_x_3d, id_dim_xu_3d, id_dim_y_xy, id_dim_yv_xy, &
249                    id_dim_y_xz, id_dim_yv_xz, id_dim_y_yz, id_dim_yv_yz, &
250                    id_dim_y_3d, id_dim_yv_3d, id_dim_zs_xy, id_dim_zs_xz, &
251                    id_dim_zs_yz, id_dim_zs_3d, id_dim_zu_xy, id_dim_zu1_xy, &
252                    id_dim_zu_xz, id_dim_zu_yz, id_dim_zu_3d, id_dim_zw_xy, &
253                    id_dim_zw_xz, id_dim_zw_yz, id_dim_zw_3d, id_set_xy, &
254                    id_set_xz, id_set_yz, id_set_3d, id_var_ind_x_yz, &
255                    id_var_ind_y_xz, id_var_ind_z_xy, id_var_time_xy, &
256                    id_var_time_xz, id_var_time_yz, id_var_time_3d, id_var_x_xy, &
257                    id_var_xu_xy, id_var_x_xz, id_var_xu_xz, id_var_x_yz, &
258                    id_var_xu_yz, id_var_x_3d, id_var_xu_3d, id_var_y_xy, &
259                    id_var_yv_xy, id_var_y_xz, id_var_yv_xz, id_var_y_yz, &
260                    id_var_yv_yz, id_var_y_3d, id_var_yv_3d, id_var_zs_xy, &
261                    id_var_zs_xz, id_var_zs_yz, id_var_zs_3d, id_var_zusi_xy, &
262                    id_var_zusi_3d, id_var_zu_xy, id_var_zu1_xy, id_var_zu_xz, &
263                    id_var_zu_yz, id_var_zu_3d, id_var_zwwi_xy, id_var_zwwi_3d, &
264                    id_var_zw_xy, id_var_zw_xz, id_var_zw_yz, id_var_zw_3d
265
266    INTEGER ::  netcdf_data_format = 2  !< NetCDF3 64bit offset format
267    INTEGER ::  netcdf_deflate = 0      !< NetCDF compression, default: no
268                                        !< compression
269
270    INTEGER(iwp), DIMENSION(10)  ::  id_var_dospx, id_var_dospy
271    INTEGER(iwp), DIMENSION(20)  ::  id_var_prt
272    INTEGER(iwp), DIMENSION(11)  ::  nc_precision
273    INTEGER(iwp), DIMENSION(dopr_norm_num) ::  id_var_norm_dopr
274
275    INTEGER(iwp), DIMENSION(dopts_num,0:10) ::  id_var_dopts
276    INTEGER(iwp), DIMENSION(0:1,100)        ::  id_var_do2d, id_var_do3d
277    INTEGER(iwp), DIMENSION(100,0:9)        ::  id_dim_z_pr, id_var_dopr, &
278                                                id_var_z_pr
279    INTEGER(iwp), DIMENSION(dots_max,0:9)   ::  id_var_dots
280
281!
282!-- Masked output
283    CHARACTER (LEN=7), DIMENSION(max_masks,0:1,100) ::  domask_unit
284
285    LOGICAL ::  output_for_t0 = .FALSE.
286
287    INTEGER(iwp), DIMENSION(1:max_masks,0:1) ::  id_dim_time_mask, id_dim_x_mask, &
288                   id_dim_xu_mask, id_dim_y_mask, id_dim_yv_mask, id_dim_zs_mask, &
289                   id_dim_zu_mask, id_dim_zw_mask, &
290                   id_set_mask, &
291                   id_var_time_mask, id_var_x_mask, id_var_xu_mask, &
292                   id_var_y_mask, id_var_yv_mask, id_var_zs_mask, &
293                   id_var_zu_mask, id_var_zw_mask, &
294                   id_var_zusi_mask, id_var_zwwi_mask
295
296    INTEGER(iwp), DIMENSION(1:max_masks,0:1,100) ::  id_var_domask
297
298
299    PUBLIC  domask_unit, dopr_unit, dopts_num, dots_label, dots_max, dots_num, &
300            dots_rad, dots_soil, dots_unit, do2d_unit, do3d_unit, id_set_mask, &
301            id_set_pr, id_set_prt, id_set_pts, id_set_sp, id_set_ts,           &
302            id_set_xy, id_set_xz, id_set_yz, id_set_3d, id_var_domask,         &
303            id_var_dopr, id_var_dopts, id_var_dospx, id_var_dospy,             &
304            id_var_dots, id_var_do2d, id_var_do3d, id_var_norm_dopr,           &
[1786]305            id_var_time_mask, id_var_time_pr, id_var_time_pts, id_var_time_sp, &
306            id_var_time_ts, id_var_time_xy, id_var_time_xz, id_var_time_yz,    &
307            id_var_time_3d, nc_stat, netcdf_data_format,                       &
308            netcdf_data_format_string, netcdf_deflate, netcdf_precision,       &
309            output_for_t0
[1783]310
311    SAVE
312
313    INTERFACE netcdf_create_dim
314       MODULE PROCEDURE netcdf_create_dim
315    END INTERFACE netcdf_create_dim
316
317    INTERFACE netcdf_create_file
318       MODULE PROCEDURE netcdf_create_file
319    END INTERFACE netcdf_create_file
320
321    INTERFACE netcdf_create_var
322       MODULE PROCEDURE netcdf_create_var
323    END INTERFACE netcdf_create_var
324
325    INTERFACE netcdf_define_header
326       MODULE PROCEDURE netcdf_define_header
327    END INTERFACE netcdf_define_header
328
329    INTERFACE netcdf_handle_error
330       MODULE PROCEDURE netcdf_handle_error
331    END INTERFACE netcdf_handle_error
332
333    INTERFACE netcdf_open_write_file
334       MODULE PROCEDURE netcdf_open_write_file
335    END INTERFACE netcdf_open_write_file
336
337    PUBLIC netcdf_create_file, netcdf_define_header, netcdf_handle_error,      &
338           netcdf_open_write_file
339
340 CONTAINS
341
342 SUBROUTINE netcdf_define_header( callmode, extend, av )
[1682]343 
[1]344#if defined( __netcdf )
345
[1691]346    USE arrays_3d,                                                             &
[1320]347        ONLY:  zu, zw
348
[1691]349    USE constants,                                                             &
[1320]350        ONLY:  pi
351
[1691]352    USE control_parameters,                                                    &
353        ONLY:  averaging_interval, averaging_interval_pr,                      &
[1833]354               data_output_pr,  domask,  dopr_n,        &
[1691]355               dopr_time_count, dopts_time_count, dots_time_count,             &
[1833]356               do2d, do2d_xz_time_count, do3d,                &
[1745]357               do2d_yz_time_count, dt_data_output_av, dt_do2d_xy, dt_do2d_xz,  &
358               dt_do2d_yz, dt_do3d, mask_size, do2d_xy_time_count,             &
359               do3d_time_count, domask_time_count, end_time, mask_i_global,    &
[1783]360               mask_j_global, mask_k_global, message_string, mid, ntdim_2d_xy, &
[1691]361               ntdim_2d_xz, ntdim_2d_yz, ntdim_3d, nz_do3d, prt_time_count,    &
[1745]362               run_description_header, section, simulated_time,                &
363               simulated_time_at_begin, skip_time_data_output_av,              &
364               skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,        &
365               skip_time_do3d, topography
[1320]366
[1691]367    USE grid_variables,                                                        &
[1320]368        ONLY:  dx, dy, zu_s_inner, zw_w_inner
369
[1691]370    USE indices,                                                               &
[1320]371        ONLY:  nx, ny, nz ,nzb, nzt
372
373    USE kinds
374
[1551]375    USE land_surface_model_mod,                                                &
[1783]376        ONLY: land_surface, nzb_soil, nzt_soil, nzs, zs
[1551]377
[1]378    USE pegrid
379
[1691]380    USE particle_attributes,                                                   &
[1320]381        ONLY:  maximum_number_of_particles, number_of_particle_groups
[1]382
[1691]383    USE profil_parameter,                                                      &
384        ONLY:  crmax, cross_profiles, dopr_index, profile_columns, profile_rows
[1320]385
[1833]386    USE spectra_mod,                                                           &
387        ONLY:  averaging_interval_sp, comp_spectra_level, data_output_sp, dosp_time_count, spectra_direction
[1320]388
[1691]389    USE statistics,                                                            &
[1320]390        ONLY:  hom, statistic_regions
391
392
[1]393    IMPLICIT NONE
394
[1682]395    CHARACTER (LEN=2)              ::  suffix                !<
396    CHARACTER (LEN=2), INTENT (IN) ::  callmode              !<
397    CHARACTER (LEN=3)              ::  suffix1               !<
398    CHARACTER (LEN=4)              ::  grid_x                !<
399    CHARACTER (LEN=4)              ::  grid_y                !<
400    CHARACTER (LEN=4)              ::  grid_z                !<
401    CHARACTER (LEN=6)              ::  mode                  !<
402    CHARACTER (LEN=10)             ::  netcdf_var_name       !<
403    CHARACTER (LEN=10)             ::  precision             !<
404    CHARACTER (LEN=10)             ::  var                   !<
405    CHARACTER (LEN=80)             ::  time_average_text     !<
406    CHARACTER (LEN=2000)           ::  char_cross_profiles   !<
407    CHARACTER (LEN=2000)           ::  var_list              !<
408    CHARACTER (LEN=2000)           ::  var_list_old          !<
[1]409
[1682]410    CHARACTER (LEN=100), DIMENSION(1:crmax) ::  cross_profiles_adj   !<
411    CHARACTER (LEN=100), DIMENSION(1:crmax) ::  cross_profiles_char  !<
[1]412
[1682]413    INTEGER(iwp) ::  av                                      !<
414    INTEGER(iwp) ::  cross_profiles_count                    !<
415    INTEGER(iwp) ::  cross_profiles_maxi                     !<
416    INTEGER(iwp) ::  delim                                   !<
417    INTEGER(iwp) ::  delim_old                               !<
418    INTEGER(iwp) ::  file_id                                 !<
419    INTEGER(iwp) ::  i                                       !<
420    INTEGER(iwp) ::  id_last                                 !<
421    INTEGER(iwp) ::  id_x                                    !<
422    INTEGER(iwp) ::  id_y                                    !<
423    INTEGER(iwp) ::  id_z                                    !<
424    INTEGER(iwp) ::  j                                       !<
425    INTEGER(iwp) ::  k                                       !<
426    INTEGER(iwp) ::  kk                                      !<
427    INTEGER(iwp) ::  ns                                      !<
428    INTEGER(iwp) ::  ns_do                                   !< actual value of ns for soil model data
429    INTEGER(iwp) ::  ns_old                                  !<
[1745]430    INTEGER(iwp) ::  ntime_count                             !< number of time levels found in file
[1682]431    INTEGER(iwp) ::  nz_old                                  !<
[951]432
[1682]433    INTEGER(iwp), SAVE ::  oldmode                           !<
[1308]434
[1682]435    INTEGER(iwp), DIMENSION(1) ::  id_dim_time_old           !<
436    INTEGER(iwp), DIMENSION(1) ::  id_dim_x_yz_old           !<
437    INTEGER(iwp), DIMENSION(1) ::  id_dim_y_xz_old           !<
438    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_sp_old          !<
439    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_xy_old          !<
440    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_3d_old          !<
441    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_mask_old        !<
[1]442
443
[1682]444    INTEGER(iwp), DIMENSION(1:crmax) ::  cross_profiles_numb !<
[951]445
[1682]446    LOGICAL ::  found                                        !<
[1]447
[1682]448    LOGICAL, INTENT (INOUT) ::  extend                       !<
[1]449
[1682]450    LOGICAL, SAVE ::  init_netcdf = .FALSE.                  !<
[1]451
[1745]452    REAL(wp), DIMENSION(1) ::  last_time_coordinate          !< last time value in file
[1]453
[1682]454    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  netcdf_data    !<
455    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  netcdf_data_2d !<
[1320]456
[1]457!
[1783]458!-- Initializing actions
[1]459    IF ( .NOT. init_netcdf )  THEN
460!
[1031]461!--    Check and set accuracy for netCDF output. First set default value
[1]462       nc_precision = NF90_REAL4
463
464       i = 1
465       DO  WHILE ( netcdf_precision(i) /= ' ' )
466          j = INDEX( netcdf_precision(i), '_' )
467          IF ( j == 0 )  THEN
[274]468             WRITE ( message_string, * ) 'netcdf_precision must contain a ', &
469                                         '"_"netcdf_precision(', i, ')="',   &
[257]470                                         TRIM( netcdf_precision(i) ),'"'
[1783]471             CALL message( 'netcdf_define_header', 'PA0241', 2, 2, 0, 6, 0 )
[1]472          ENDIF
473
474          var       = netcdf_precision(i)(1:j-1)
475          precision = netcdf_precision(i)(j+1:)
476
477          IF ( precision == 'NF90_REAL4' )  THEN
478             j = NF90_REAL4
479          ELSEIF ( precision == 'NF90_REAL8' )  THEN
480             j = NF90_REAL8
481          ELSE
[257]482             WRITE ( message_string, * ) 'illegal netcdf precision: ',  &
483                                         'netcdf_precision(', i, ')="', &
484                                         TRIM( netcdf_precision(i) ),'"'
[1783]485             CALL message( 'netcdf_define_header', 'PA0242', 1, 2, 0, 6, 0 )
[1]486          ENDIF
487
488          SELECT CASE ( var )
489             CASE ( 'xy' )
490                nc_precision(1) = j
491             CASE ( 'xz' )
492                nc_precision(2) = j
493             CASE ( 'yz' )
494                nc_precision(3) = j
495             CASE ( '2d' )
496                nc_precision(1:3) = j
497             CASE ( '3d' )
498                nc_precision(4) = j
499             CASE ( 'pr' )
500                nc_precision(5) = j
501             CASE ( 'ts' )
502                nc_precision(6) = j
503             CASE ( 'sp' )
504                nc_precision(7) = j
505             CASE ( 'prt' )
506                nc_precision(8) = j
[410]507             CASE ( 'masks' )
508                nc_precision(11) = j
[1]509             CASE ( 'all' )
510                nc_precision    = j
511
512             CASE DEFAULT
[274]513                WRITE ( message_string, * ) 'unknown variable in inipar ',    & 
514                                  'assignment: netcdf_precision(', i, ')="',  &
[257]515                                            TRIM( netcdf_precision(i) ),'"'
[1783]516                CALL message( 'netcdf_define_header', 'PA0243', 1, 2, 0, 6, 0 )
[1]517
518          END SELECT
519
520          i = i + 1
[410]521          IF ( i > 50 )  EXIT
[1]522       ENDDO
523
[1783]524!
525!--    Check for allowed parameter range
526       IF ( netcdf_deflate < 0  .OR.  netcdf_deflate > 9 )  THEN
527          WRITE ( message_string, '(A,I3,A)' ) 'netcdf_deflate out of ' //     &
528                                      'range & given value: ', netcdf_deflate, &
529                                      ', allowed range: 0-9'
530          CALL message( 'netcdf_define_header', 'PA0355', 2, 2, 0, 6, 0 )
531       ENDIF
532!
533!--    Data compression does not work with parallel NetCDF/HDF5
534       IF ( netcdf_deflate > 0  .AND.  netcdf_data_format /= 3 )  THEN
535          message_string = 'netcdf_deflate reset to 0'
536          CALL message( 'netcdf_define_header', 'PA0356', 0, 1, 0, 6, 0 )
537
538          netcdf_deflate = 0
539       ENDIF
540
[1]541       init_netcdf = .TRUE.
542
543    ENDIF
544
545!
546!-- Determine the mode to be processed
547    IF ( extend )  THEN
548       mode = callmode // '_ext'
549    ELSE
550       mode = callmode // '_new'
551    ENDIF
552
553!
[410]554!-- Select the mode to be processed. Possibilities are 3d, mask, xy, xz, yz,
[292]555!-- pr and ts.
[1]556    SELECT CASE ( mode )
557
[410]558       CASE ( 'ma_new' )
559
560!
561!--       decompose actual parameter file_id (=formal parameter av) into
562!--       mid and av
563          file_id = av
[564]564          IF ( file_id <= 200+max_masks )  THEN
565             mid = file_id - 200
[410]566             av = 0
567          ELSE
[564]568             mid = file_id - (200+max_masks)
[410]569             av = 1
570          ENDIF
571
572!
573!--       Define some global attributes of the dataset
574          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
[1783]575                                  'Conventions', 'COARDS' )
576          CALL netcdf_handle_error( 'netcdf_define_header', 464 )
[410]577
578          IF ( av == 0 )  THEN
579             time_average_text = ' '
580          ELSE
581             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
582                                                            averaging_interval
583          ENDIF
584          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'title', &
585                                  TRIM( run_description_header ) //    &
586                                  TRIM( time_average_text ) )
[1783]587          CALL netcdf_handle_error( 'netcdf_define_header', 465 )
[410]588          IF ( av == 1 )  THEN
589             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
590             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
591                                     'time_avg', TRIM( time_average_text ) )
[1783]592             CALL netcdf_handle_error( 'netcdf_define_header', 466 )
[410]593          ENDIF
594
595!
596!--       Define time coordinate for volume data (unlimited dimension)
[1783]597          CALL netcdf_create_dim( id_set_mask(mid,av), 'time', NF90_UNLIMITED, &
598                                  id_dim_time_mask(mid,av), 467 )
599          CALL netcdf_create_var( id_set_mask(mid,av),                         &
600                                  (/ id_dim_time_mask(mid,av) /), 'time',      &
601                                  NF90_DOUBLE, id_var_time_mask(mid,av),       &
602                                 'seconds', '', 468, 469, 000 )
[410]603!
604!--       Define spatial dimensions and coordinates:
605!--       Define vertical coordinate grid (zu grid)
[1783]606          CALL netcdf_create_dim( id_set_mask(mid,av), 'zu_3d',                &
607                                  mask_size(mid,3), id_dim_zu_mask(mid,av),    &
608                                  470 )
609          CALL netcdf_create_var( id_set_mask(mid,av),                         &
610                                  (/ id_dim_zu_mask(mid,av) /), 'zu_3d',       &
611                                  NF90_DOUBLE, id_var_zu_mask(mid,av),         &
612                                  'meters', '', 471, 472, 000 )
[410]613!
614!--       Define vertical coordinate grid (zw grid)
[1783]615          CALL netcdf_create_dim( id_set_mask(mid,av), 'zw_3d',                &
616                                  mask_size(mid,3), id_dim_zw_mask(mid,av),    &
617                                  473 )
618          CALL netcdf_create_var( id_set_mask(mid,av),                         &
619                                  (/ id_dim_zw_mask(mid,av) /), 'zw_3d',       &
620                                  NF90_DOUBLE, id_var_zw_mask(mid,av),         &
621                                 'meters', '', 474, 475, 000 )
[410]622!
623!--       Define x-axis (for scalar position)
[1783]624          CALL netcdf_create_dim( id_set_mask(mid,av), 'x', mask_size(mid,1),  &
625                                  id_dim_x_mask(mid,av), 476 )
626          CALL netcdf_create_var( id_set_mask(mid,av),                         &
627                                  (/ id_dim_x_mask(mid,av) /), 'x',            &
628                                  NF90_DOUBLE, id_var_x_mask(mid,av),          &
629                                  'meters', '', 477, 478, 000 )
[410]630!
631!--       Define x-axis (for u position)
[1783]632          CALL netcdf_create_dim( id_set_mask(mid,av), 'xu', mask_size(mid,1), &
633                                  id_dim_xu_mask(mid,av), 479 )
634          CALL netcdf_create_var( id_set_mask(mid,av),                         &
635                                  (/ id_dim_xu_mask(mid,av) /), 'xu',          &
636                                  NF90_DOUBLE, id_var_xu_mask(mid,av),         &
637                                  'meters', '', 480, 481, 000 )
[410]638!
639!--       Define y-axis (for scalar position)
[1783]640          CALL netcdf_create_dim( id_set_mask(mid,av), 'y', mask_size(mid,2),  &
641                                  id_dim_y_mask(mid,av), 482 )
642          CALL netcdf_create_var( id_set_mask(mid,av),                         &
643                                  (/ id_dim_y_mask(mid,av) /), 'y',            &
644                                  NF90_DOUBLE, id_var_y_mask(mid,av),          &
645                                  'meters', '', 483, 484, 000 )
[410]646!
647!--       Define y-axis (for v position)
[1783]648          CALL netcdf_create_dim( id_set_mask(mid,av), 'yv', mask_size(mid,2), &
649                                  id_dim_yv_mask(mid,av), 485 )
650          CALL netcdf_create_var( id_set_mask(mid,av),                         &
651                                  (/ id_dim_yv_mask(mid,av) /),                &
652                                  'yv', NF90_DOUBLE, id_var_yv_mask(mid,av),   &
653                                  'meters', '', 486, 487, 000 )
[410]654!
655!--       In case of non-flat topography define 2d-arrays containing the height
[1551]656!--       information
[410]657          IF ( TRIM( topography ) /= 'flat' )  THEN
658!
659!--          Define zusi = zu(nzb_s_inner)
[1783]660             CALL netcdf_create_var( id_set_mask(mid,av),                      &
661                                     (/ id_dim_x_mask(mid,av),                 &
662                                        id_dim_y_mask(mid,av) /), 'zusi',      &
663                                     NF90_DOUBLE, id_var_zusi_mask(mid,av),    &
664                                     'meters', 'zu(nzb_s_inner)', 488, 489,    &
665                                     490 )
[410]666!             
667!--          Define zwwi = zw(nzb_w_inner)
[1783]668             CALL netcdf_create_var( id_set_mask(mid,av),                      &
669                                     (/ id_dim_x_mask(mid,av),                 &
670                                        id_dim_y_mask(mid,av) /), 'zwwi',      &
671                                     NF90_DOUBLE, id_var_zwwi_mask(mid,av),    &
672                                     'meters', 'zw(nzb_w_inner)', 491, 492,    &
673                                     493 )
[410]674          ENDIF             
[1551]675 
676          IF ( land_surface )  THEN
677!
678!--          Define vertical coordinate grid (zw grid)
[1783]679             CALL netcdf_create_dim( id_set_mask(mid,av), 'zs_3d',             &
680                                     mask_size(mid,3), id_dim_zs_mask(mid,av), &
681                                     536 )
682             CALL netcdf_create_var( id_set_mask(mid,av),                      &
683                                     (/ id_dim_zs_mask(mid,av) /), 'zs_3d',    &
684                                     NF90_DOUBLE, id_var_zs_mask(mid,av),      &
685                                     'meters', '', 537, 555, 000 )
[1551]686          ENDIF
687
[410]688!
689!--       Define the variables
690          var_list = ';'
691          i = 1
692
693          DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
694
695!
696!--          Check for the grid
697             found = .TRUE.
698             SELECT CASE ( domask(mid,av,i) )
699!
700!--             Most variables are defined on the scalar grid
[1691]701                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
702                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',        &
703                       'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr',            &
704                       'rad_sw_hr', 'rho', 's', 'sa', 'vpt'  )
[410]705
706                   grid_x = 'x'
707                   grid_y = 'y'
708                   grid_z = 'zu'
709!
710!--             u grid
711                CASE ( 'u' )
712
713                   grid_x = 'xu'
714                   grid_y = 'y'
715                   grid_z = 'zu'
716!
717!--             v grid
718                CASE ( 'v' )
719
720                   grid_x = 'x'
721                   grid_y = 'yv'
722                   grid_z = 'zu'
723!
724!--             w grid
[1691]725                CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',   &
726                       'w' )
[410]727
728                   grid_x = 'x'
729                   grid_y = 'y'
730                   grid_z = 'zw'
[1551]731!
732!--             soil grid
[1691]733                CASE ( 'm_soil', 't_soil' )
[410]734
[1551]735                   grid_x = 'x'
736                   grid_y = 'y'
737                   grid_z = 'zs'
738
[410]739                CASE DEFAULT
740!
741!--                Check for user-defined quantities
[1783]742                   CALL user_define_netcdf_grid( domask(mid,av,i), found,      &
[410]743                                                 grid_x, grid_y, grid_z )
744
[1783]745                   IF ( .NOT. found )  THEN
746                      WRITE ( message_string, * ) 'no grid defined for',       &
747                           ' variable ', TRIM( domask(mid,av,i) )
748                      CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, &
749                                    6, 0 )
750                   ENDIF
751
[410]752             END SELECT
753
754!
755!--          Select the respective dimension ids
756             IF ( grid_x == 'x' )  THEN
757                id_x = id_dim_x_mask(mid,av)
758             ELSEIF ( grid_x == 'xu' )  THEN
759                id_x = id_dim_xu_mask(mid,av)
760             ENDIF
761
762             IF ( grid_y == 'y' )  THEN
763                id_y = id_dim_y_mask(mid,av)
764             ELSEIF ( grid_y == 'yv' )  THEN
765                id_y = id_dim_yv_mask(mid,av)
766             ENDIF
767
768             IF ( grid_z == 'zu' )  THEN
769                id_z = id_dim_zu_mask(mid,av)
770             ELSEIF ( grid_z == 'zw' )  THEN
771                id_z = id_dim_zw_mask(mid,av)
[1551]772             ELSEIF ( grid_z == "zs" )  THEN
773                id_z = id_dim_zs_mask(mid,av)
[410]774             ENDIF
775
776!
777!--          Define the grid
[1783]778             CALL netcdf_create_var( id_set_mask(mid,av), (/ id_x, id_y, id_z, &
779                                     id_dim_time_mask(mid,av) /),              &
780                                     domask(mid,av,i), nc_precision(11),       &
781                                     id_var_domask(mid,av,i),                  &
782                                     TRIM( domask_unit(mid,av,i) ),            &
783                                     domask(mid,av,i), 494, 495, 496 )
[410]784
785             var_list = TRIM( var_list ) // TRIM( domask(mid,av,i) ) // ';'
786
787             i = i + 1
788
789          ENDDO
790
791!
792!--       No arrays to output
793          IF ( i == 1 )  RETURN
794
795!
796!--       Write the list of variables as global attribute (this is used by
797!--       restart runs and by combine_plot_fields)
798          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
799                                  'VAR_LIST', var_list )
[1783]800          CALL netcdf_handle_error( 'netcdf_define_header', 497 )
[410]801
802!
[1031]803!--       Leave netCDF define mode
[410]804          nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
[1783]805          CALL netcdf_handle_error( 'netcdf_define_header', 498 )
[410]806
807!
808!--       Write data for x (shifted by +dx/2) and xu axis
809          ALLOCATE( netcdf_data(mask_size(mid,1)) )
810
[1353]811          netcdf_data = ( mask_i_global(mid,:mask_size(mid,1)) + 0.5_wp ) * dx
[410]812
813          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_x_mask(mid,av), &
814                                  netcdf_data, start = (/ 1 /),               &
815                                  count = (/ mask_size(mid,1) /) )
[1783]816          CALL netcdf_handle_error( 'netcdf_define_header', 499 )
[410]817
818          netcdf_data = mask_i_global(mid,:mask_size(mid,1)) * dx
819
820          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_xu_mask(mid,av),&
821                                  netcdf_data, start = (/ 1 /),               &
822                                  count = (/ mask_size(mid,1) /) )
[1783]823          CALL netcdf_handle_error( 'netcdf_define_header', 500 )
[410]824
825          DEALLOCATE( netcdf_data )
826
827!
828!--       Write data for y (shifted by +dy/2) and yv axis
829          ALLOCATE( netcdf_data(mask_size(mid,2)) )
830
[1353]831          netcdf_data = ( mask_j_global(mid,:mask_size(mid,2)) + 0.5_wp ) * dy
[410]832
833          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_y_mask(mid,av), &
834                                  netcdf_data, start = (/ 1 /),               &
835                                  count = (/ mask_size(mid,2) /))
[1783]836          CALL netcdf_handle_error( 'netcdf_define_header', 501 )
[410]837
838          netcdf_data = mask_j_global(mid,:mask_size(mid,2)) * dy
839
840          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_yv_mask(mid,av), &
841                                  netcdf_data, start = (/ 1 /),    &
842                                  count = (/ mask_size(mid,2) /))
[1783]843          CALL netcdf_handle_error( 'netcdf_define_header', 502 )
[410]844
845          DEALLOCATE( netcdf_data )
846
847!
848!--       Write zu and zw data (vertical axes)
849          ALLOCATE( netcdf_data(mask_size(mid,3)) )
850
851          netcdf_data = zu( mask_k_global(mid,:mask_size(mid,3)) )
852
853          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
854                                  netcdf_data, start = (/ 1 /), &
855                                  count = (/ mask_size(mid,3) /) )
[1783]856          CALL netcdf_handle_error( 'netcdf_define_header', 503 )
[410]857
858          netcdf_data = zw( mask_k_global(mid,:mask_size(mid,3)) )
859
860          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
861                                  netcdf_data, start = (/ 1 /), &
862                                  count = (/ mask_size(mid,3) /) )
[1783]863          CALL netcdf_handle_error( 'netcdf_define_header', 504 )
[410]864
865          DEALLOCATE( netcdf_data )
866
867!
868!--       In case of non-flat topography write height information
869          IF ( TRIM( topography ) /= 'flat' )  THEN
870
871             ALLOCATE( netcdf_data_2d(mask_size(mid,1),mask_size(mid,2)) )
872             netcdf_data_2d = zu_s_inner( mask_i_global(mid,:mask_size(mid,1)),&
873                                          mask_j_global(mid,:mask_size(mid,2)) )
874
875             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),         &
876                                     id_var_zusi_mask(mid,av),    &
877                                     netcdf_data_2d,              &
878                                     start = (/ 1, 1 /),          &
879                                     count = (/ mask_size(mid,1), &
880                                                mask_size(mid,2) /) )
[1783]881             CALL netcdf_handle_error( 'netcdf_define_header', 505 )
[410]882
883             netcdf_data_2d = zw_w_inner( mask_i_global(mid,:mask_size(mid,1)),&
884                                          mask_j_global(mid,:mask_size(mid,2)) )
885
886             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),         &
887                                     id_var_zwwi_mask(mid,av),    &
888                                     netcdf_data_2d,              &
889                                     start = (/ 1, 1 /),          &
890                                     count = (/ mask_size(mid,1), &
891                                                mask_size(mid,2) /) )
[1783]892             CALL netcdf_handle_error( 'netcdf_define_header', 506 )
[410]893
894             DEALLOCATE( netcdf_data_2d )
895
896          ENDIF
[1551]897
898          IF ( land_surface )  THEN
[410]899!
[1551]900!--          Write zs data (vertical axes for soil model), use negative values
901!--          to indicate soil depth
902             ALLOCATE( netcdf_data(mask_size(mid,3)) )
903
904             netcdf_data = zs( mask_k_global(mid,:mask_size(mid,3)) )
905
906             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zs_mask(mid,av), &
907                                     netcdf_data, start = (/ 1 /), &
908                                     count = (/ mask_size(mid,3) /) )
[1783]909             CALL netcdf_handle_error( 'netcdf_define_header', 538 )
[1551]910
911             DEALLOCATE( netcdf_data )
912
913          ENDIF
914
915!
[410]916!--       restore original parameter file_id (=formal parameter av) into av
917          av = file_id
918
919
920       CASE ( 'ma_ext' )
921
922!
923!--       decompose actual parameter file_id (=formal parameter av) into
924!--       mid and av
925          file_id = av
[564]926          IF ( file_id <= 200+max_masks )  THEN
927             mid = file_id - 200
[410]928             av = 0
929          ELSE
[564]930             mid = file_id - (200+max_masks)
[410]931             av = 1
932          ENDIF
933
934!
935!--       Get the list of variables and compare with the actual run.
936!--       First var_list_old has to be reset, since GET_ATT does not assign
937!--       trailing blanks.
938          var_list_old = ' '
939          nc_stat = NF90_GET_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'VAR_LIST',&
940                                  var_list_old )
[1783]941          CALL netcdf_handle_error( 'netcdf_define_header', 507 )
[410]942
943          var_list = ';'
944          i = 1
945          DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
946             var_list = TRIM(var_list) // TRIM( domask(mid,av,i) ) // ';'
947             i = i + 1
948          ENDDO
949
950          IF ( av == 0 )  THEN
951             var = '(mask)'
952          ELSE
953             var = '(mask_av)'
954          ENDIF
955
956          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
[1031]957             WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ), &
[410]958                  ' data for mask', mid, ' from previous run found,', &
959                  '&but this file cannot be extended due to variable ', &
960                  'mismatch.&New file is created instead.'
[564]961             CALL message( 'define_netcdf_header', 'PA0335', 0, 1, 0, 6, 0 )
[410]962             extend = .FALSE.
963             RETURN
964          ENDIF
965
966!
967!--       Get and compare the number of vertical gridpoints
[1596]968          nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'zu_3d', &
[410]969                                    id_var_zu_mask(mid,av) )
[1783]970          CALL netcdf_handle_error( 'netcdf_define_header', 508 )
[410]971
972          nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av),     &
973                                           id_var_zu_mask(mid,av),  &
974                                           dimids = id_dim_zu_mask_old )
[1783]975          CALL netcdf_handle_error( 'netcdf_define_header', 509 )
[410]976          id_dim_zu_mask(mid,av) = id_dim_zu_mask_old(1)
977
978          nc_stat = NF90_INQUIRE_DIMENSION( id_set_mask(mid,av),    &
979                                            id_dim_zu_mask(mid,av), &
980                                            len = nz_old )
[1783]981          CALL netcdf_handle_error( 'netcdf_define_header', 510 )
[410]982
983          IF ( mask_size(mid,3) /= nz_old )  THEN
[1031]984             WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ), &
[410]985                  ' data for mask', mid, ' from previous run found,', &
986                  '&but this file cannot be extended due to mismatch in ', &
987                  ' number of&vertical grid points.', &
988                  '&New file is created instead.'
[564]989             CALL message( 'define_netcdf_header', 'PA0336', 0, 1, 0, 6, 0 )
[410]990             extend = .FALSE.
991             RETURN
992          ENDIF
993
994!
995!--       Get the id of the time coordinate (unlimited coordinate) and its
996!--       last index on the file. The next time level is plmask..count+1.
997!--       The current time must be larger than the last output time
998!--       on the file.
999          nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'time', &
1000                                    id_var_time_mask(mid,av) )
[1783]1001          CALL netcdf_handle_error( 'netcdf_define_header', 511 )
[410]1002
1003          nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av), &
1004                                           id_var_time_mask(mid,av), &
1005                                           dimids = id_dim_time_old )
[1783]1006          CALL netcdf_handle_error( 'netcdf_define_header', 512 )
[410]1007          id_dim_time_mask(mid,av) = id_dim_time_old(1)
1008
1009          nc_stat = NF90_INQUIRE_DIMENSION( id_set_mask(mid,av), &
1010                                            id_dim_time_mask(mid,av), &
1011                                            len = domask_time_count(mid,av) )
[1783]1012          CALL netcdf_handle_error( 'netcdf_define_header', 513 )
[410]1013
1014          nc_stat = NF90_GET_VAR( id_set_mask(mid,av), &
1015                                  id_var_time_mask(mid,av), &
1016                                  last_time_coordinate,              &
1017                                  start = (/ domask_time_count(mid,av) /), &
1018                                  count = (/ 1 /) )
[1783]1019          CALL netcdf_handle_error( 'netcdf_define_header', 514 )
[410]1020
1021          IF ( last_time_coordinate(1) >= simulated_time )  THEN
[1031]1022             WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ), &
[410]1023                  ' data for mask', mid, ' from previous run found,', &
1024                  '&but this file cannot be extended because the current ', &
1025                  'output time&is less or equal than the last output time ', &
1026                  'on this file.&New file is created instead.'
[564]1027             CALL message( 'define_netcdf_header', 'PA0337', 0, 1, 0, 6, 0 )
[410]1028             domask_time_count(mid,av) = 0
1029             extend = .FALSE.
1030             RETURN
1031          ENDIF
1032
1033!
1034!--       Dataset seems to be extendable.
1035!--       Now get the variable ids.
1036          i = 1
1037          DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
1038             nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), &
1039                                       TRIM( domask(mid,av,i) ), &
1040                                       id_var_domask(mid,av,i) )
[1783]1041             CALL netcdf_handle_error( 'netcdf_define_header', 515 )
[410]1042             i = i + 1
1043          ENDDO
1044
1045!
1046!--       Update the title attribute on file
1047!--       In order to avoid 'data mode' errors if updated attributes are larger
1048!--       than their original size, NF90_PUT_ATT is called in 'define mode'
1049!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
1050!--       performance loss due to data copying; an alternative strategy would be
1051!--       to ensure equal attribute size in a job chain. Maybe revise later.
1052          IF ( av == 0 )  THEN
1053             time_average_text = ' '
1054          ELSE
1055             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
1056                                                            averaging_interval
1057          ENDIF
1058          nc_stat = NF90_REDEF( id_set_mask(mid,av) )
[1783]1059          CALL netcdf_handle_error( 'netcdf_define_header', 516 )
[410]1060          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'title', &
1061                                  TRIM( run_description_header ) //    &
1062                                  TRIM( time_average_text ) )
[1783]1063          CALL netcdf_handle_error( 'netcdf_define_header', 517 )
[410]1064          nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
[1783]1065          CALL netcdf_handle_error( 'netcdf_define_header', 518 )
[1031]1066          WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ), &
[410]1067               ' data for mask', mid, ' from previous run found.', &
1068               '&This file will be extended.'
[564]1069          CALL message( 'define_netcdf_header', 'PA0338', 0, 0, 0, 6, 0 )
[410]1070!
1071!--       restore original parameter file_id (=formal parameter av) into av
1072          av = file_id
1073
1074
[1]1075       CASE ( '3d_new' )
1076
1077!
1078!--       Define some global attributes of the dataset
1079          nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'Conventions', &
1080                                  'COARDS' )
[1783]1081          CALL netcdf_handle_error( 'netcdf_define_header', 62 )
[1]1082
1083          IF ( av == 0 )  THEN
1084             time_average_text = ' '
1085          ELSE
1086             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
1087                                                            averaging_interval
1088          ENDIF
1089          nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'title', &
1090                                  TRIM( run_description_header ) //    &
1091                                  TRIM( time_average_text ) )
[1783]1092          CALL netcdf_handle_error( 'netcdf_define_header', 63 )
[1]1093          IF ( av == 1 )  THEN
1094             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
1095             nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'time_avg', &
1096                                     TRIM( time_average_text ) )
[1783]1097             CALL netcdf_handle_error( 'netcdf_define_header', 63 )
[1]1098          ENDIF
1099
1100!
[1308]1101!--       Define time coordinate for volume data.
1102!--       For parallel output the time dimensions has to be limited, otherwise
1103!--       the performance drops significantly.
1104          IF ( netcdf_data_format < 5 )  THEN
[1783]1105             CALL netcdf_create_dim( id_set_3d(av), 'time', NF90_UNLIMITED,    &
1106                                     id_dim_time_3d(av), 64 )
[1308]1107          ELSE
[1783]1108             CALL netcdf_create_dim( id_set_3d(av), 'time', ntdim_3d(av),      &
1109                                     id_dim_time_3d(av), 523 )
[1308]1110          ENDIF
[1]1111
[1783]1112          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_time_3d(av) /),     &
1113                                  'time', NF90_DOUBLE, id_var_time_3d(av),     &
1114                                  'seconds', '', 65, 66, 00 )
[1]1115!
1116!--       Define spatial dimensions and coordinates:
1117!--       Define vertical coordinate grid (zu grid)
[1783]1118          CALL netcdf_create_dim( id_set_3d(av), 'zu_3d', nz_do3d-nzb+1,       &
1119                                  id_dim_zu_3d(av), 67 )
1120          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zu_3d(av) /),       &
1121                                  'zu_3d', NF90_DOUBLE, id_var_zu_3d(av),      &
1122                                  'meters', '', 68, 69, 00 )
[1]1123!
1124!--       Define vertical coordinate grid (zw grid)
[1783]1125          CALL netcdf_create_dim( id_set_3d(av), 'zw_3d', nz_do3d-nzb+1,       &
1126                                  id_dim_zw_3d(av), 70 )
1127          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zw_3d(av) /),       &
1128                                  'zw_3d', NF90_DOUBLE, id_var_zw_3d(av),      &
1129                                  'meters', '', 71, 72, 00 )
[1]1130!
1131!--       Define x-axis (for scalar position)
[1783]1132          CALL netcdf_create_dim( id_set_3d(av), 'x', nx+2, id_dim_x_3d(av),   &
1133                                  73 )
1134          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av) /), 'x',   &
1135                                  NF90_DOUBLE, id_var_x_3d(av), 'meters', '',  &
1136                                  74, 75, 00 )
[1]1137!
1138!--       Define x-axis (for u position)
[1783]1139          CALL netcdf_create_dim( id_set_3d(av), 'xu', nx+2, id_dim_xu_3d(av), &
1140                                  358 )
1141          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_xu_3d(av) /), 'xu', &
1142                                  NF90_DOUBLE, id_var_xu_3d(av), 'meters', '', &
1143                                  359, 360, 000 )
[1]1144!
1145!--       Define y-axis (for scalar position)
[1783]1146          CALL netcdf_create_dim( id_set_3d(av), 'y', ny+2, id_dim_y_3d(av),   &
1147                                  76 )
1148          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_y_3d(av) /), 'y',   &
1149                                  NF90_DOUBLE, id_var_y_3d(av), 'meters', '',  &
1150                                  77, 78, 00 )
[1]1151!
1152!--       Define y-axis (for v position)
[1783]1153          CALL netcdf_create_dim( id_set_3d(av), 'yv', ny+2, id_dim_yv_3d(av), &
1154                                  361 )
1155          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_yv_3d(av) /), 'yv', &
1156                                  NF90_DOUBLE, id_var_yv_3d(av), 'meters', '', &
1157                                  362, 363, 000 )
[1]1158!
[48]1159!--       In case of non-flat topography define 2d-arrays containing the height
[1551]1160!--       information
[48]1161          IF ( TRIM( topography ) /= 'flat' )  THEN
1162!
1163!--          Define zusi = zu(nzb_s_inner)
[1783]1164             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av),        &
1165                                     id_dim_y_3d(av) /), 'zusi', NF90_DOUBLE,  &
1166                                     id_var_zusi_3d(av), 'meters',             &
1167                                     'zu(nzb_s_inner)', 413, 414, 415 )
[48]1168!             
1169!--          Define zwwi = zw(nzb_w_inner)
[1783]1170             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av),        &
1171                                     id_dim_y_3d(av) /), 'zwwi', NF90_DOUBLE,  &
1172                                     id_var_zwwi_3d(av), 'meters',             &
1173                                     'zw(nzb_w_inner)', 416, 417, 418 )
[48]1174
1175          ENDIF             
1176
[1551]1177          IF ( land_surface )  THEN
1178!
1179!--          Define vertical coordinate grid (zs grid)
[1783]1180             CALL netcdf_create_dim( id_set_3d(av), 'zs_3d',                   &
1181                                     nzt_soil-nzb_soil+1, id_dim_zs_3d(av), 70 )
1182             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zs_3d(av) /),    &
1183                                     'zs_3d', NF90_DOUBLE, id_var_zs_3d(av),   &
1184                                     'meters', '', 71, 72, 00 )
[48]1185
[1551]1186          ENDIF
1187
[48]1188!
[1]1189!--       Define the variables
1190          var_list = ';'
1191          i = 1
1192
1193          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
1194
1195!
1196!--          Check for the grid
1197             found = .TRUE.
1198             SELECT CASE ( do3d(av,i) )
1199!
1200!--             Most variables are defined on the scalar grid
[1691]1201                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
1202                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho', &
1203                       's', 'sa', 'vpt' , 'rad_lw_cs_hr', 'rad_lw_hr',         &
1204                       'rad_sw_cs_hr', 'rad_sw_hr' )
[1]1205
1206                   grid_x = 'x'
1207                   grid_y = 'y'
1208                   grid_z = 'zu'
1209!
1210!--             u grid
1211                CASE ( 'u' )
1212
1213                   grid_x = 'xu'
1214                   grid_y = 'y'
1215                   grid_z = 'zu'
1216!
1217!--             v grid
1218                CASE ( 'v' )
1219
1220                   grid_x = 'x'
1221                   grid_y = 'yv'
1222                   grid_z = 'zu'
1223!
1224!--             w grid
[1691]1225                CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',   &
1226                       'w' )
[1]1227
1228                   grid_x = 'x'
1229                   grid_y = 'y'
1230                   grid_z = 'zw'
[1551]1231!
1232!--             soil grid
[1691]1233                CASE ( 'm_soil', 't_soil' )
[1]1234
[1551]1235                   grid_x = 'x'
1236                   grid_y = 'y'
1237                   grid_z = 'zs'
1238
[1]1239                CASE DEFAULT
1240!
1241!--                Check for user-defined quantities
1242                   CALL user_define_netcdf_grid( do3d(av,i), found, grid_x, &
1243                                                 grid_y, grid_z )
1244
[1783]1245                   IF ( .NOT. found )  THEN
1246                      WRITE ( message_string, * ) 'no grid defined for varia', &
1247                                                  'ble ', TRIM( do3d(av,i) )
1248                      CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, &
1249                                    6, 0 )
1250                   ENDIF
1251
[1]1252             END SELECT
1253
1254!
1255!--          Select the respective dimension ids
1256             IF ( grid_x == 'x' )  THEN
1257                id_x = id_dim_x_3d(av)
1258             ELSEIF ( grid_x == 'xu' )  THEN
1259                id_x = id_dim_xu_3d(av)
1260             ENDIF
1261
1262             IF ( grid_y == 'y' )  THEN
1263                id_y = id_dim_y_3d(av)
1264             ELSEIF ( grid_y == 'yv' )  THEN
1265                id_y = id_dim_yv_3d(av)
1266             ENDIF
1267
1268             IF ( grid_z == 'zu' )  THEN
1269                id_z = id_dim_zu_3d(av)
1270             ELSEIF ( grid_z == 'zw' )  THEN
1271                id_z = id_dim_zw_3d(av)
[1551]1272             ELSEIF ( grid_z == 'zs' )  THEN
1273                id_z = id_dim_zs_3d(av)
[1]1274             ENDIF
1275
1276!
1277!--          Define the grid
[1783]1278             CALL netcdf_create_var( id_set_3d(av),(/ id_x, id_y, id_z,        &
1279                                     id_dim_time_3d(av) /), do3d(av,i),        &
1280                                     nc_precision(4), id_var_do3d(av,i),       &
1281                                     TRIM( do3d_unit(av,i) ), do3d(av,i), 79,  &
1282                                     80, 357 )
[1031]1283#if defined( __netcdf4_parallel )
[1308]1284             IF ( netcdf_data_format > 4 )  THEN
[493]1285!
[1308]1286!--             Set no fill for every variable to increase performance.
1287                nc_stat = NF90_DEF_VAR_FILL( id_set_3d(av),     &
1288                                             id_var_do3d(av,i), &
1289                                             1, 0 )
[1783]1290                CALL netcdf_handle_error( 'netcdf_define_header', 532 )
[1308]1291!
1292!--             Set collective io operations for parallel io
[493]1293                nc_stat = NF90_VAR_PAR_ACCESS( id_set_3d(av),     &
1294                                               id_var_do3d(av,i), &
1295                                               NF90_COLLECTIVE )
[1783]1296                CALL netcdf_handle_error( 'netcdf_define_header', 445 )
[493]1297             ENDIF
1298#endif
[1783]1299             var_list = TRIM( var_list ) // TRIM( do3d(av,i) ) // ';'
1300
[1]1301             i = i + 1
1302
1303          ENDDO
1304
1305!
1306!--       No arrays to output
1307          IF ( i == 1 )  RETURN
1308
1309!
1310!--       Write the list of variables as global attribute (this is used by
1311!--       restart runs and by combine_plot_fields)
1312          nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
1313                                  var_list )
[1783]1314          CALL netcdf_handle_error( 'netcdf_define_header', 81 )
[1]1315
1316!
[1308]1317!--       Set general no fill, otherwise the performance drops significantly for
1318!--       parallel output.
1319          nc_stat = NF90_SET_FILL( id_set_3d(av), NF90_NOFILL, oldmode )
[1783]1320          CALL netcdf_handle_error( 'netcdf_define_header', 528 )
[1308]1321
1322!
[1031]1323!--       Leave netCDF define mode
[1]1324          nc_stat = NF90_ENDDEF( id_set_3d(av) )
[1783]1325          CALL netcdf_handle_error( 'netcdf_define_header', 82 )
[1]1326
1327!
[1308]1328!--       These data are only written by PE0 for parallel output to increase
1329!--       the performance.
1330          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
1331!
1332!--          Write data for x (shifted by +dx/2) and xu axis
1333             ALLOCATE( netcdf_data(0:nx+1) )
[1]1334
[1308]1335             DO  i = 0, nx+1
1336                netcdf_data(i) = ( i + 0.5 ) * dx
1337             ENDDO
[1]1338
[1308]1339             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_x_3d(av),  &
1340                                     netcdf_data, start = (/ 1 /),    &
1341                                     count = (/ nx+2 /) )
[1783]1342             CALL netcdf_handle_error( 'netcdf_define_header', 83 )
[1]1343
[1308]1344             DO  i = 0, nx+1
1345                netcdf_data(i) = i * dx
1346             ENDDO
[1]1347
[1308]1348             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_xu_3d(av), &
1349                                     netcdf_data, start = (/ 1 /),    &
1350                                     count = (/ nx+2 /) )
[1783]1351             CALL netcdf_handle_error( 'netcdf_define_header', 385 )
[1]1352
[1308]1353             DEALLOCATE( netcdf_data )
[1]1354
1355!
[1308]1356!--          Write data for y (shifted by +dy/2) and yv axis
1357             ALLOCATE( netcdf_data(0:ny+1) )
[1]1358
[1308]1359             DO  i = 0, ny+1
[1353]1360                netcdf_data(i) = ( i + 0.5_wp ) * dy
[1308]1361             ENDDO
[1]1362
[1308]1363             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_y_3d(av),  &
1364                                     netcdf_data, start = (/ 1 /),    &
1365                                     count = (/ ny+2 /) )
[1783]1366             CALL netcdf_handle_error( 'netcdf_define_header', 84 )
[1]1367
[1308]1368             DO  i = 0, ny+1
1369                netcdf_data(i) = i * dy
1370             ENDDO
[1]1371
[1308]1372             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_yv_3d(av), &
1373                                     netcdf_data, start = (/ 1 /),    &
1374                                     count = (/ ny+2 /))
[1783]1375             CALL netcdf_handle_error( 'netcdf_define_header', 387 )
[1]1376
[1308]1377             DEALLOCATE( netcdf_data )
[1]1378
1379!
[1308]1380!--          Write zu and zw data (vertical axes)
1381             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zu_3d(av),  &
1382                                     zu(nzb:nz_do3d), start = (/ 1 /), &
1383                                     count = (/ nz_do3d-nzb+1 /) )
[1783]1384             CALL netcdf_handle_error( 'netcdf_define_header', 85 )
[1]1385
[263]1386
[1308]1387             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zw_3d(av),  &
1388                                     zw(nzb:nz_do3d), start = (/ 1 /), &
1389                                     count = (/ nz_do3d-nzb+1 /) )
[1783]1390             CALL netcdf_handle_error( 'netcdf_define_header', 86 )
[1]1391
[48]1392!
[1308]1393!--          In case of non-flat topography write height information
1394             IF ( TRIM( topography ) /= 'flat' )  THEN
[48]1395
[1308]1396                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av), &
1397                                        zu_s_inner(0:nx+1,0:ny+1), &
1398                                        start = (/ 1, 1 /), &
1399                                        count = (/ nx+2, ny+2 /) )
[1783]1400                CALL netcdf_handle_error( 'netcdf_define_header', 419 )
[48]1401
[1308]1402                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av), &
1403                                        zw_w_inner(0:nx+1,0:ny+1), &
1404                                        start = (/ 1, 1 /), &
1405                                        count = (/ nx+2, ny+2 /) )
[1783]1406                CALL netcdf_handle_error( 'netcdf_define_header', 420 )
[48]1407
[1308]1408             ENDIF
[1551]1409
1410             IF ( land_surface )  THEN
1411!
1412!--             Write zs grid
1413                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zs_3d(av),  &
1414                                        - zs(nzb_soil:nzt_soil), start = (/ 1 /), &
1415                                        count = (/ nzt_soil-nzb_soil+1 /) )
[1783]1416                CALL netcdf_handle_error( 'netcdf_define_header', 86 )
[1551]1417             ENDIF
1418
[48]1419          ENDIF
1420
[1]1421       CASE ( '3d_ext' )
1422
1423!
1424!--       Get the list of variables and compare with the actual run.
1425!--       First var_list_old has to be reset, since GET_ATT does not assign
1426!--       trailing blanks.
1427          var_list_old = ' '
1428          nc_stat = NF90_GET_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
1429                                  var_list_old )
[1783]1430          CALL netcdf_handle_error( 'netcdf_define_header', 87 )
[1]1431
1432          var_list = ';'
1433          i = 1
1434          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
1435             var_list = TRIM(var_list) // TRIM( do3d(av,i) ) // ';'
1436             i = i + 1
1437          ENDDO
1438
1439          IF ( av == 0 )  THEN
1440             var = '(3d)'
1441          ELSE
1442             var = '(3d_av)'
1443          ENDIF
1444
1445          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
[1031]1446             message_string = 'netCDF file for volume data ' //             &
[292]1447                              TRIM( var ) // ' from previous run found,' // &
[257]1448                              '&but this file cannot be extended due to' // &
1449                              ' variable mismatch.' //                      &
1450                              '&New file is created instead.'
1451             CALL message( 'define_netcdf_header', 'PA0245', 0, 1, 0, 6, 0 )
[1]1452             extend = .FALSE.
1453             RETURN
1454          ENDIF
1455
1456!
1457!--       Get and compare the number of vertical gridpoints
1458          nc_stat = NF90_INQ_VARID( id_set_3d(av), 'zu_3d', id_var_zu_3d(av) )
[1783]1459          CALL netcdf_handle_error( 'netcdf_define_header', 88 )
[1]1460
1461          nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_zu_3d(av), &
1462                                           dimids = id_dim_zu_3d_old )
[1783]1463          CALL netcdf_handle_error( 'netcdf_define_header', 89 )
[1]1464          id_dim_zu_3d(av) = id_dim_zu_3d_old(1)
1465
1466          nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_zu_3d(av), &
1467                                            len = nz_old )
[1783]1468          CALL netcdf_handle_error( 'netcdf_define_header', 90 )
[1]1469
1470          IF ( nz_do3d-nzb+1 /= nz_old )  THEN
[1031]1471              message_string = 'netCDF file for volume data ' //             &
[292]1472                               TRIM( var ) // ' from previous run found,' // &
[257]1473                               '&but this file cannot be extended due to' // &
1474                               ' mismatch in number of' //                   &
1475                               '&vertical grid points (nz_do3d).' //         &
1476                               '&New file is created instead.'
1477             CALL message( 'define_netcdf_header', 'PA0246', 0, 1, 0, 6, 0 )
[1]1478             extend = .FALSE.
1479             RETURN
1480          ENDIF
1481
1482!
1483!--       Get the id of the time coordinate (unlimited coordinate) and its
1484!--       last index on the file. The next time level is pl3d..count+1.
1485!--       The current time must be larger than the last output time
1486!--       on the file.
1487          nc_stat = NF90_INQ_VARID( id_set_3d(av), 'time', id_var_time_3d(av) )
[1783]1488          CALL netcdf_handle_error( 'netcdf_define_header', 91 )
[1]1489
1490          nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_time_3d(av), &
1491                                           dimids = id_dim_time_old )
[1783]1492          CALL netcdf_handle_error( 'netcdf_define_header', 92 )
[263]1493
[1]1494          id_dim_time_3d(av) = id_dim_time_old(1)
1495
1496          nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_time_3d(av), &
[1308]1497                                            len = ntime_count )
[1783]1498          CALL netcdf_handle_error( 'netcdf_define_header', 93 )
[1]1499
[1308]1500!
1501!--       For non-parallel output use the last output time level of the netcdf
1502!--       file because the time dimension is unlimited. In case of parallel
1503!--       output the variable ntime_count could get the value of 9*10E36 because
1504!--       the time dimension is limited.
1505          IF ( netcdf_data_format < 5 ) do3d_time_count(av) = ntime_count
1506
[1]1507          nc_stat = NF90_GET_VAR( id_set_3d(av), id_var_time_3d(av), &
1508                                  last_time_coordinate,              &
1509                                  start = (/ do3d_time_count(av) /), &
1510                                  count = (/ 1 /) )
[1783]1511          CALL netcdf_handle_error( 'netcdf_define_header', 94 )
[1]1512
1513          IF ( last_time_coordinate(1) >= simulated_time )  THEN
[1031]1514             message_string = 'netCDF file for volume data ' // &
[292]1515                              TRIM( var ) // ' from previous run found,' // &
[257]1516                              '&but this file cannot be extended becaus' // &
1517                              'e the current output time' //                &
1518                              '&is less or equal than the last output t' // &
1519                              'ime on this file.' //                        &
1520                              '&New file is created instead.'
1521             CALL message( 'define_netcdf_header', 'PA0247', 0, 1, 0, 6, 0 )
[1]1522             do3d_time_count(av) = 0
1523             extend = .FALSE.
1524             RETURN
1525          ENDIF
1526
[1308]1527          IF ( netcdf_data_format > 4 )  THEN
[1745]1528!
1529!--          Set needed time levels (ntdim) to
1530!--          saved time levels + to be saved time levels.
1531             IF ( av == 0 )  THEN
1532                ntdim_3d(0) = do3d_time_count(0) + CEILING(                &
1533                              ( end_time - MAX( skip_time_do3d,            &
1534                                                simulated_time_at_begin )  &
1535                              ) / dt_do3d )
1536             ELSE
1537                ntdim_3d(1) = do3d_time_count(1) + CEILING(                &
1538                              ( end_time - MAX( skip_time_data_output_av,  &
1539                                                simulated_time_at_begin )  &
1540                              ) / dt_data_output_av )
1541             ENDIF
1542
1543!
1544!--          Check if the needed number of output time levels is increased
1545!--          compared to the number of time levels in the existing file.
[1308]1546             IF ( ntdim_3d(av) > ntime_count )  THEN
1547                message_string = 'netCDF file for volume data ' // &
1548                                 TRIM( var ) // ' from previous run found,' // &
1549                                 '&but this file cannot be extended becaus' // &
1550                                 'e the number of output time levels&has b' // &
1551                                 'een increased compared to the previous s' // &
[1745]1552                                 'imulation.' //                               &
[1308]1553                                 '&New file is created instead.'
1554                CALL message( 'define_netcdf_header', 'PA0388', 0, 1, 0, 6, 0 )
1555                do3d_time_count(av) = 0
1556                extend = .FALSE.
[1745]1557!
1558!--             Recalculate the needed time levels for the new file.
1559                IF ( av == 0 )  THEN
1560                   ntdim_3d(0) = CEILING(                               &
1561                           ( end_time - MAX( skip_time_do3d,            &
1562                                             simulated_time_at_begin )  &
1563                           ) / dt_do3d )
1564                ELSE
1565                   ntdim_3d(1) = CEILING(                               &
1566                           ( end_time - MAX( skip_time_data_output_av,  &
1567                                             simulated_time_at_begin )  &
1568                           ) / dt_data_output_av )
1569                ENDIF
[1308]1570                RETURN
1571             ENDIF
1572          ENDIF
1573
[1]1574!
1575!--       Dataset seems to be extendable.
1576!--       Now get the variable ids.
1577          i = 1
1578          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
1579             nc_stat = NF90_INQ_VARID( id_set_3d(av), TRIM( do3d(av,i) ), &
1580                                       id_var_do3d(av,i) )
[1783]1581             CALL netcdf_handle_error( 'netcdf_define_header', 95 )
[1031]1582#if defined( __netcdf4_parallel )
[493]1583!
1584!--          Set collective io operations for parallel io
[1031]1585             IF ( netcdf_data_format > 4 )  THEN
[493]1586                nc_stat = NF90_VAR_PAR_ACCESS( id_set_3d(av),     &
1587                                               id_var_do3d(av,i), &
1588                                               NF90_COLLECTIVE )
[1783]1589                CALL netcdf_handle_error( 'netcdf_define_header', 453 )
[493]1590             ENDIF
1591#endif
[1]1592             i = i + 1
1593          ENDDO
1594
1595!
[359]1596!--       Update the title attribute on file
1597!--       In order to avoid 'data mode' errors if updated attributes are larger
1598!--       than their original size, NF90_PUT_ATT is called in 'define mode'
1599!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
1600!--       performance loss due to data copying; an alternative strategy would be
1601!--       to ensure equal attribute size. Maybe revise later.
1602          IF ( av == 0 )  THEN
1603             time_average_text = ' '
1604          ELSE
1605             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
1606                                                            averaging_interval
1607          ENDIF
1608          nc_stat = NF90_REDEF( id_set_3d(av) )
[1783]1609          CALL netcdf_handle_error( 'netcdf_define_header', 429 )
[1]1610          nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'title', &
[359]1611                                  TRIM( run_description_header ) //    &
1612                                  TRIM( time_average_text ) )
[1783]1613          CALL netcdf_handle_error( 'netcdf_define_header', 96 )
[359]1614          nc_stat = NF90_ENDDEF( id_set_3d(av) )
[1783]1615          CALL netcdf_handle_error( 'netcdf_define_header', 430 )
[1031]1616          message_string = 'netCDF file for volume data ' //             &
[257]1617                           TRIM( var ) // ' from previous run found.' // &
1618                           '&This file will be extended.'
1619          CALL message( 'define_netcdf_header', 'PA0248', 0, 0, 0, 6, 0 )
[1]1620
1621       CASE ( 'xy_new' )
1622
1623!
1624!--       Define some global attributes of the dataset
1625          nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'Conventions', &
1626                                  'COARDS' )
[1783]1627          CALL netcdf_handle_error( 'netcdf_define_header', 97 )
[1]1628
1629          IF ( av == 0 )  THEN
1630             time_average_text = ' '
1631          ELSE
1632             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
1633                                                            averaging_interval
1634          ENDIF
1635          nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'title', &
1636                                  TRIM( run_description_header ) //    &
1637                                  TRIM( time_average_text ) )
[1783]1638          CALL netcdf_handle_error( 'netcdf_define_header', 98 )
[1]1639          IF ( av == 1 )  THEN
1640             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
1641             nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'time_avg', &
1642                                     TRIM( time_average_text ) )
[1783]1643             CALL netcdf_handle_error( 'netcdf_define_header', 98 )
[1]1644          ENDIF
1645
1646!
[1308]1647!--       Define time coordinate for xy sections.
1648!--       For parallel output the time dimensions has to be limited, otherwise
1649!--       the performance drops significantly.
1650          IF ( netcdf_data_format < 5 )  THEN
[1783]1651             CALL netcdf_create_dim( id_set_xy(av), 'time', NF90_UNLIMITED,    &
1652                                     id_dim_time_xy(av), 99 )
[1308]1653          ELSE
[1783]1654             CALL netcdf_create_dim( id_set_xy(av), 'time', ntdim_2d_xy(av),   &
1655                                     id_dim_time_xy(av), 524 )
[1308]1656          ENDIF
[1]1657
[1783]1658          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_time_xy(av) /),     &
1659                                  'time', NF90_DOUBLE, id_var_time_xy(av),     &
1660                                  'seconds', '', 100, 101, 000 )
[1]1661!
1662!--       Define the spatial dimensions and coordinates for xy-sections.
1663!--       First, determine the number of horizontal sections to be written.
1664          IF ( section(1,1) == -9999 )  THEN
1665             RETURN
1666          ELSE
1667             ns = 1
1668             DO WHILE ( section(ns,1) /= -9999  .AND.  ns <= 100 )
1669                ns = ns + 1
1670             ENDDO
1671             ns = ns - 1
1672          ENDIF
1673
1674!
1675!--       Define vertical coordinate grid (zu grid)
[1783]1676          CALL netcdf_create_dim( id_set_xy(av), 'zu_xy', ns,                  &
1677                                  id_dim_zu_xy(av), 102 )
1678          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu_xy(av) /),       &
1679                                  'zu_xy', NF90_DOUBLE, id_var_zu_xy(av),      &
1680                                  'meters', '', 103, 104, 000 )
[1]1681!
1682!--       Define vertical coordinate grid (zw grid)
[1783]1683          CALL netcdf_create_dim( id_set_xy(av), 'zw_xy', ns,                  &
1684                                  id_dim_zw_xy(av), 105 )
1685          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zw_xy(av) /),       &
1686                                  'zw_xy', NF90_DOUBLE, id_var_zw_xy(av),      &
1687                                  'meters', '', 106, 107, 000 )
[1]1688
[1551]1689          IF ( land_surface )  THEN
1690
1691             ns_do = 0
1692             DO WHILE ( section(ns_do+1,1) < nzs )
1693                ns_do = ns_do + 1
1694             ENDDO
[1]1695!
[1551]1696!--          Define vertical coordinate grid (zs grid)
[1783]1697             CALL netcdf_create_dim( id_set_xy(av), 'zs_xy', ns_do,            &
1698                                     id_dim_zs_xy(av), 539 )
1699             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zs_xy(av) /),    &
1700                                     'zs_xy', NF90_DOUBLE, id_var_zs_xy(av),   &
1701                                     'meters', '', 540, 541, 000 )
[1551]1702
1703          ENDIF
1704
1705!
[1]1706!--       Define a pseudo vertical coordinate grid for the surface variables
1707!--       u* and t* to store their height level
[1783]1708          CALL netcdf_create_dim( id_set_xy(av), 'zu1_xy', 1,                  &
1709                                  id_dim_zu1_xy(av), 108 )
1710          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu1_xy(av) /),      &
1711                                  'zu1_xy', NF90_DOUBLE, id_var_zu1_xy(av),    &
1712                                  'meters', '', 109, 110, 000 )
[1]1713!
1714!--       Define a variable to store the layer indices of the horizontal cross
1715!--       sections, too
[1783]1716          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu_xy(av) /),       &
1717                                  'ind_z_xy', NF90_DOUBLE,                     &
1718                                  id_var_ind_z_xy(av), 'gridpoints', '', 111,  &
1719                                  112, 000 )
[1]1720!
1721!--       Define x-axis (for scalar position)
[1783]1722          CALL netcdf_create_dim( id_set_xy(av), 'x', nx+2, id_dim_x_xy(av),   &
1723                                  113 )
1724          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av) /), 'x',   &
1725                                  NF90_DOUBLE, id_var_x_xy(av), 'meters', '',  &
1726                                  114, 115, 000 )
[1]1727!
1728!--       Define x-axis (for u position)
[1783]1729          CALL netcdf_create_dim( id_set_xy(av), 'xu', nx+2,                   &
1730                                  id_dim_xu_xy(av), 388 )
1731          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_xu_xy(av) /), 'xu', &
1732                                  NF90_DOUBLE, id_var_xu_xy(av), 'meters', '', &
1733                                  389, 390, 000 )
[1]1734!
1735!--       Define y-axis (for scalar position)
[1783]1736          CALL netcdf_create_dim( id_set_xy(av), 'y', ny+2, id_dim_y_xy(av),   &
1737                                  116 )
1738          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_y_xy(av) /), 'y',   &
1739                                  NF90_DOUBLE, id_var_y_xy(av), 'meters', '',  &
1740                                  117, 118, 000 )
[1]1741!
1742!--       Define y-axis (for scalar position)
[1783]1743          CALL netcdf_create_dim( id_set_xy(av), 'yv', ny+2,                   &
1744                                  id_dim_yv_xy(av), 364 )
1745          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_yv_xy(av) /), 'yv', &
1746                                  NF90_DOUBLE, id_var_yv_xy(av), 'meters', '', &
1747                                  365, 366, 000 )
[1]1748!
[48]1749!--       In case of non-flat topography define 2d-arrays containing the height
[1551]1750!--       information
[48]1751          IF ( TRIM( topography ) /= 'flat' )  THEN
1752!
1753!--          Define zusi = zu(nzb_s_inner)
[1783]1754             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),        &
1755                                     id_dim_y_xy(av) /), 'zusi', NF90_DOUBLE,  &
1756                                     id_var_zusi_xy(av), 'meters',             &
1757                                     'zu(nzb_s_inner)', 421, 422, 423 )
[48]1758!             
1759!--          Define zwwi = zw(nzb_w_inner)
[1783]1760             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),        &
1761                                     id_dim_y_xy(av) /), 'zwwi', NF90_DOUBLE,  &
1762                                     id_var_zwwi_xy(av), 'meters',             &
1763                                     'zw(nzb_w_inner)', 424, 425, 426 )
[48]1764
1765          ENDIF
1766
1767!
[1]1768!--       Define the variables
1769          var_list = ';'
1770          i = 1
1771
1772          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
1773
1774             IF ( INDEX( do2d(av,i), 'xy' ) /= 0 )  THEN
1775!
1776!--             If there is a star in the variable name (u* or t*), it is a
1777!--             surface variable. Define it with id_dim_zu1_xy.
1778                IF ( INDEX( do2d(av,i), '*' ) /= 0 )  THEN
1779
[1783]1780                   CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),  &
1781                                           id_dim_y_xy(av), id_dim_zu1_xy(av), &
1782                                           id_dim_time_xy(av) /), do2d(av,i),  &
1783                                           nc_precision(1), id_var_do2d(av,i), &
1784                                           TRIM( do2d_unit(av,i) ),            &
1785                                           do2d(av,i), 119, 120, 354 )
[1]1786
1787                ELSE
1788
1789!
1790!--                Check for the grid
1791                   found = .TRUE.
1792                   SELECT CASE ( do2d(av,i) )
1793!
1794!--                   Most variables are defined on the zu grid
[1691]1795                      CASE ( 'e_xy', 'lpt_xy', 'nr_xy', 'p_xy', 'pc_xy',       &
1796                             'pr_xy', 'prr_xy', 'pt_xy', 'q_xy', 'qc_xy',      &
1797                             'ql_xy', 'ql_c_xy', 'ql_v_xy', 'ql_vp_xy',        &
1798                             'qr_xy', 'qv_xy', 'rad_lw_cs_hr_xy',              &
1799                             'rad_lw_hr_xy', 'rad_sw_cs_hr_xy', 'rad_sw_hr_xy',& 
1800                             'rho_xy', 's_xy', 'sa_xy', 'vpt_xy' )
[1]1801
1802                         grid_x = 'x'
1803                         grid_y = 'y'
1804                         grid_z = 'zu'
1805!
1806!--                   u grid
1807                      CASE ( 'u_xy' )
1808
1809                         grid_x = 'xu'
1810                         grid_y = 'y'
1811                         grid_z = 'zu'
1812!
1813!--                   v grid
1814                      CASE ( 'v_xy' )
1815
1816                         grid_x = 'x'
1817                         grid_y = 'yv'
1818                         grid_z = 'zu'
1819!
1820!--                   w grid
[1691]1821                      CASE ( 'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy',  &
1822                             'rad_sw_out_xy' , 'w_xy' )
[1]1823
1824                         grid_x = 'x'
1825                         grid_y = 'y'
1826                         grid_z = 'zw'
[1551]1827!
1828!--                   soil grid
[1691]1829                      CASE ( 'm_soil_xy', 't_soil_xy' )
[1551]1830                         grid_x = 'x'
1831                         grid_y = 'y'
1832                         grid_z = 'zs'
[1]1833
1834                      CASE DEFAULT
1835!
1836!--                      Check for user-defined quantities
1837                         CALL user_define_netcdf_grid( do2d(av,i), found, &
1838                                                       grid_x, grid_y, grid_z )
1839
[1783]1840                         IF ( .NOT. found )  THEN
1841                            WRITE ( message_string, * ) 'no grid defined for', &
1842                                                ' variable ', TRIM( do2d(av,i) )
1843                            CALL message( 'define_netcdf_header', 'PA0244',    &
1844                                          0, 1, 0, 6, 0 )
1845                         ENDIF
1846
[1]1847                   END SELECT
1848
1849!
1850!--                Select the respective dimension ids
1851                   IF ( grid_x == 'x' )  THEN
1852                      id_x = id_dim_x_xy(av)
1853                   ELSEIF ( grid_x == 'xu' )  THEN
1854                      id_x = id_dim_xu_xy(av)
1855                   ENDIF
1856
1857                   IF ( grid_y == 'y' )  THEN
1858                      id_y = id_dim_y_xy(av)
1859                   ELSEIF ( grid_y == 'yv' )  THEN
1860                      id_y = id_dim_yv_xy(av)
1861                   ENDIF
1862
1863                   IF ( grid_z == 'zu' )  THEN
1864                      id_z = id_dim_zu_xy(av)
1865                   ELSEIF ( grid_z == 'zw' )  THEN
1866                      id_z = id_dim_zw_xy(av)
[1551]1867                   ELSEIF ( grid_z == 'zs' )  THEN
1868                      id_z = id_dim_zs_xy(av)
[1]1869                   ENDIF
1870
1871!
1872!--                Define the grid
[1783]1873                   CALL netcdf_create_var( id_set_xy(av), (/ id_x, id_y, id_z, &
1874                                           id_dim_time_xy(av) /), do2d(av,i),  &
1875                                           nc_precision(1), id_var_do2d(av,i), &
1876                                           TRIM( do2d_unit(av,i) ),            &
1877                                           do2d(av,i), 119, 120, 354 )
[1]1878
1879                ENDIF
1880
[1031]1881#if defined( __netcdf4_parallel )
[1308]1882                IF ( netcdf_data_format > 4 )  THEN
[493]1883!
[1308]1884!--                Set no fill for every variable to increase performance.
1885                   nc_stat = NF90_DEF_VAR_FILL( id_set_xy(av),     &
1886                                                id_var_do2d(av,i), &
1887                                                1, 0 )
[1783]1888                   CALL netcdf_handle_error( 'netcdf_define_header', 533 )
[1308]1889!
1890!--                Set collective io operations for parallel io
[493]1891                   nc_stat = NF90_VAR_PAR_ACCESS( id_set_xy(av),     &
1892                                                  id_var_do2d(av,i), &
1893                                                  NF90_COLLECTIVE )
[1783]1894                   CALL netcdf_handle_error( 'netcdf_define_header', 448 )
[493]1895                ENDIF
1896#endif
[1783]1897                var_list = TRIM( var_list) // TRIM( do2d(av,i) ) // ';'
1898
[1]1899             ENDIF
1900
1901             i = i + 1
1902
1903          ENDDO
1904
1905!
1906!--       No arrays to output. Close the netcdf file and return.
1907          IF ( i == 1 )  RETURN
1908
1909!
1910!--       Write the list of variables as global attribute (this is used by
1911!--       restart runs and by combine_plot_fields)
1912          nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'VAR_LIST', &
1913                                  var_list )
[1783]1914          CALL netcdf_handle_error( 'netcdf_define_header', 121 )
[1]1915
1916!
[1308]1917!--       Set general no fill, otherwise the performance drops significantly for
1918!--       parallel output.
1919          nc_stat = NF90_SET_FILL( id_set_xy(av), NF90_NOFILL, oldmode )
[1783]1920          CALL netcdf_handle_error( 'netcdf_define_header', 529 )
[1308]1921
1922!
[1031]1923!--       Leave netCDF define mode
[1]1924          nc_stat = NF90_ENDDEF( id_set_xy(av) )
[1783]1925          CALL netcdf_handle_error( 'netcdf_define_header', 122 )
[1]1926
1927!
[1308]1928!--       These data are only written by PE0 for parallel output to increase
1929!--       the performance.
1930          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
[1]1931
1932!
[1308]1933!--          Write axis data: z_xy, x, y
1934             ALLOCATE( netcdf_data(1:ns) )
[1]1935
1936!
[1308]1937!--          Write zu data
1938             DO  i = 1, ns
1939                IF( section(i,1) == -1 )  THEN
[1353]1940                   netcdf_data(i) = -1.0_wp  ! section averaged along z
[1308]1941                ELSE
1942                   netcdf_data(i) = zu( section(i,1) )
1943                ENDIF
1944             ENDDO
1945             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu_xy(av), &
1946                                     netcdf_data, start = (/ 1 /),    &
1947                                     count = (/ ns /) )
[1783]1948             CALL netcdf_handle_error( 'netcdf_define_header', 123 )
[1]1949
1950!
[1308]1951!--          Write zw data
1952             DO  i = 1, ns
1953                IF( section(i,1) == -1 )  THEN
[1353]1954                   netcdf_data(i) = -1.0_wp  ! section averaged along z
[1308]1955                ELSE
1956                   netcdf_data(i) = zw( section(i,1) )
1957                ENDIF
1958             ENDDO
1959             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zw_xy(av), &
1960                                     netcdf_data, start = (/ 1 /),    &
1961                                     count = (/ ns /) )
[1783]1962             CALL netcdf_handle_error( 'netcdf_define_header', 124 )
[1]1963
[1308]1964!
[1551]1965!--             Write zs data
1966             IF ( land_surface )  THEN
1967                ns_do = 0
1968                DO  i = 1, ns
1969                   IF( section(i,1) == -1 )  THEN
1970                      netcdf_data(i) = 1.0_wp  ! section averaged along z
1971                      ns_do = ns_do + 1
1972                   ELSEIF ( section(i,1) < nzs )  THEN
1973                      netcdf_data(i) = - zs( section(i,1) )
1974                      ns_do = ns_do + 1
1975                   ENDIF
1976                ENDDO
1977
1978                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zs_xy(av), &
1979                                        netcdf_data(1:ns_do), start = (/ 1 /),    &
1980                                        count = (/ ns_do /) )
[1783]1981                CALL netcdf_handle_error( 'netcdf_define_header', 124 )
[1551]1982
1983             ENDIF
1984
1985!
[1308]1986!--          Write gridpoint number data
1987             netcdf_data(1:ns) = section(1:ns,1)
1988             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_ind_z_xy(av), &
1989                                     netcdf_data, start = (/ 1 /),       &
1990                                     count = (/ ns /) )
[1783]1991             CALL netcdf_handle_error( 'netcdf_define_header', 125 )
[1]1992
[1308]1993             DEALLOCATE( netcdf_data )
1994
[1]1995!
[1308]1996!--          Write the cross section height u*, t*
1997             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu1_xy(av), &
1998                                     (/ zu(nzb+1) /), start = (/ 1 /), &
1999                                     count = (/ 1 /) )
[1783]2000             CALL netcdf_handle_error( 'netcdf_define_header', 126 )
[1]2001
2002!
[1308]2003!--          Write data for x (shifted by +dx/2) and xu axis
2004             ALLOCATE( netcdf_data(0:nx+1) )
[1]2005
[1308]2006             DO  i = 0, nx+1
[1353]2007                netcdf_data(i) = ( i + 0.5_wp ) * dx
[1308]2008             ENDDO
[1]2009
[1308]2010             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_x_xy(av), &
2011                                     netcdf_data, start = (/ 1 /),   &
2012                                     count = (/ nx+2 /) )
[1783]2013             CALL netcdf_handle_error( 'netcdf_define_header', 127 )
[1]2014
[1308]2015             DO  i = 0, nx+1
2016                netcdf_data(i) = i * dx
2017             ENDDO
[1]2018
[1308]2019             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_xu_xy(av), &
2020                                     netcdf_data, start = (/ 1 /),    &
2021                                     count = (/ nx+2 /) )
[1783]2022             CALL netcdf_handle_error( 'netcdf_define_header', 367 )
[1]2023
[1308]2024             DEALLOCATE( netcdf_data )
[1]2025
2026!
[1308]2027!--          Write data for y (shifted by +dy/2) and yv axis
2028             ALLOCATE( netcdf_data(0:ny+1) )
[1]2029
[1308]2030             DO  i = 0, ny+1
[1353]2031                netcdf_data(i) = ( i + 0.5_wp ) * dy
[1308]2032             ENDDO
[1]2033
[1308]2034             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_y_xy(av), &
2035                                     netcdf_data, start = (/ 1 /),   &
2036                                     count = (/ ny+2 /))
[1783]2037             CALL netcdf_handle_error( 'netcdf_define_header', 128 )
[1]2038
[1308]2039             DO  i = 0, ny+1
2040                netcdf_data(i) = i * dy
2041             ENDDO
[1]2042
[1308]2043             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_yv_xy(av), &
2044                                     netcdf_data, start = (/ 1 /),    &
2045                                     count = (/ ny+2 /))
[1783]2046             CALL netcdf_handle_error( 'netcdf_define_header', 368 )
[1]2047
[1308]2048             DEALLOCATE( netcdf_data )
[1]2049
[48]2050!
[1308]2051!--          In case of non-flat topography write height information
2052             IF ( TRIM( topography ) /= 'flat' )  THEN
[1]2053
[1308]2054                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av), &
2055                                        zu_s_inner(0:nx+1,0:ny+1), &
2056                                        start = (/ 1, 1 /), &
2057                                        count = (/ nx+2, ny+2 /) )
[1783]2058                CALL netcdf_handle_error( 'netcdf_define_header', 427 )
[48]2059
[1308]2060                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av), &
2061                                        zw_w_inner(0:nx+1,0:ny+1), &
2062                                        start = (/ 1, 1 /), &
2063                                        count = (/ nx+2, ny+2 /) )
[1783]2064                CALL netcdf_handle_error( 'netcdf_define_header', 428 )
[48]2065
[1308]2066             ENDIF
[1551]2067
2068
2069
[48]2070          ENDIF
2071
[1]2072       CASE ( 'xy_ext' )
2073
2074!
2075!--       Get the list of variables and compare with the actual run.
2076!--       First var_list_old has to be reset, since GET_ATT does not assign
2077!--       trailing blanks.
2078          var_list_old = ' '
2079          nc_stat = NF90_GET_ATT( id_set_xy(av), NF90_GLOBAL, 'VAR_LIST', &
2080                                  var_list_old )
[1783]2081          CALL netcdf_handle_error( 'netcdf_define_header', 129 )
[1]2082
2083          var_list = ';'
2084          i = 1
2085          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2086             IF ( INDEX( do2d(av,i), 'xy' ) /= 0 )  THEN
[519]2087                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
[1]2088             ENDIF
2089             i = i + 1
2090          ENDDO
2091
2092          IF ( av == 0 )  THEN
2093             var = '(xy)'
2094          ELSE
2095             var = '(xy_av)'
2096          ENDIF
2097
2098          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
[1031]2099             message_string = 'netCDF file for cross-sections ' //           &
[292]2100                              TRIM( var ) // ' from previous run found,' //  &
[257]2101                              '& but this file cannot be extended due to' // &
2102                              ' variable mismatch.' //                       &
2103                              '&New file is created instead.'
2104             CALL message( 'define_netcdf_header', 'PA0249', 0, 1, 0, 6, 0 )
[1]2105             extend = .FALSE.
2106             RETURN
2107          ENDIF
2108
2109!
2110!--       Calculate the number of current sections
2111          ns = 1
2112          DO WHILE ( section(ns,1) /= -9999  .AND.  ns <= 100 )
2113             ns = ns + 1
2114          ENDDO
2115          ns = ns - 1
2116
2117!
2118!--       Get and compare the number of horizontal cross sections
2119          nc_stat = NF90_INQ_VARID( id_set_xy(av), 'zu_xy', id_var_zu_xy(av) )
[1783]2120          CALL netcdf_handle_error( 'netcdf_define_header', 130 )
[1]2121
2122          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_zu_xy(av), &
2123                                           dimids = id_dim_zu_xy_old )
[1783]2124          CALL netcdf_handle_error( 'netcdf_define_header', 131 )
[1]2125          id_dim_zu_xy(av) = id_dim_zu_xy_old(1)
2126
2127          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_zu_xy(av), &
2128                                            len = ns_old )
[1783]2129          CALL netcdf_handle_error( 'netcdf_define_header', 132 )
[1]2130
2131          IF ( ns /= ns_old )  THEN
[1031]2132             message_string = 'netCDF file for cross-sections ' //          &
[292]2133                              TRIM( var ) // ' from previous run found,' // &
[257]2134                              '&but this file cannot be extended due to' // &
2135                              ' mismatch in number of' //                   &
2136                              '&cross sections.' //                         &
2137                              '&New file is created instead.'
2138             CALL message( 'define_netcdf_header', 'PA0250', 0, 1, 0, 6, 0 )
[1]2139             extend = .FALSE.
2140             RETURN
2141          ENDIF
2142
2143!
2144!--       Get and compare the heights of the cross sections
2145          ALLOCATE( netcdf_data(1:ns_old) )
2146
2147          nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_zu_xy(av), netcdf_data )
[1783]2148          CALL netcdf_handle_error( 'netcdf_define_header', 133 )
[1]2149
2150          DO  i = 1, ns
2151             IF ( section(i,1) /= -1 )  THEN
2152                IF ( zu(section(i,1)) /= netcdf_data(i) )  THEN
[1031]2153                   message_string = 'netCDF file for cross-sections ' //     &
[292]2154                               TRIM( var ) // ' from previous run found,' // &
[274]2155                               '&but this file cannot be extended' //        &
2156                               ' due to mismatch in cross' //                &
2157                               '&section levels.' //                         &
2158                               '&New file is created instead.'
2159                   CALL message( 'define_netcdf_header', 'PA0251', &
2160                                                                 0, 1, 0, 6, 0 )
[1]2161                   extend = .FALSE.
2162                   RETURN
2163                ENDIF
2164             ELSE
[1353]2165                IF ( -1.0_wp /= netcdf_data(i) )  THEN
[1031]2166                   message_string = 'netCDF file for cross-sections ' //     &
[292]2167                               TRIM( var ) // ' from previous run found,' // &
[274]2168                               '&but this file cannot be extended' //        &
2169                               ' due to mismatch in cross' //                &
2170                               '&section levels.' //                         &
2171                               '&New file is created instead.'
2172                   CALL message( 'define_netcdf_header', 'PA0251',&
2173                                                                 0, 1, 0, 6, 0 )
[1]2174                   extend = .FALSE.
2175                   RETURN
2176                ENDIF
2177             ENDIF
2178          ENDDO
2179
2180          DEALLOCATE( netcdf_data )
2181
2182!
2183!--       Get the id of the time coordinate (unlimited coordinate) and its
2184!--       last index on the file. The next time level is do2d..count+1.
2185!--       The current time must be larger than the last output time
2186!--       on the file.
2187          nc_stat = NF90_INQ_VARID( id_set_xy(av), 'time', id_var_time_xy(av) )
[1783]2188          CALL netcdf_handle_error( 'netcdf_define_header', 134 )
[1]2189
2190          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_time_xy(av), &
2191                                           dimids = id_dim_time_old )
[1783]2192          CALL netcdf_handle_error( 'netcdf_define_header', 135 )
[1]2193          id_dim_time_xy(av) = id_dim_time_old(1)
2194
2195          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_time_xy(av), &
[1308]2196                                            len = ntime_count )
[1783]2197          CALL netcdf_handle_error( 'netcdf_define_header', 136 )
[1]2198
[1308]2199!
2200!--       For non-parallel output use the last output time level of the netcdf
2201!--       file because the time dimension is unlimited. In case of parallel
2202!--       output the variable ntime_count could get the value of 9*10E36 because
2203!--       the time dimension is limited.
2204          IF ( netcdf_data_format < 5 ) do2d_xy_time_count(av) = ntime_count
2205
[1]2206          nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_time_xy(av),    &
2207                                  last_time_coordinate,                 &
2208                                  start = (/ do2d_xy_time_count(av) /), &
2209                                  count = (/ 1 /) )
[1783]2210          CALL netcdf_handle_error( 'netcdf_define_header', 137 )
[1]2211
2212          IF ( last_time_coordinate(1) >= simulated_time )  THEN
[1031]2213             message_string = 'netCDF file for cross sections ' //          &
[292]2214                              TRIM( var ) // ' from previous run found,' // &
[257]2215                              '&but this file cannot be extended becaus' // &
2216                              'e the current output time' //                &
2217                              '&is less or equal than the last output t' // &
2218                              'ime on this file.' //                        &
2219                              '&New file is created instead.'
2220             CALL message( 'define_netcdf_header', 'PA0252', 0, 1, 0, 6, 0 )
[1]2221             do2d_xy_time_count(av) = 0
2222             extend = .FALSE.
2223             RETURN
2224          ENDIF
2225
[1308]2226          IF ( netcdf_data_format > 4 )  THEN
[1745]2227!
2228!--          Set needed time levels (ntdim) to
2229!--          saved time levels + to be saved time levels.
2230             IF ( av == 0 )  THEN
2231                ntdim_2d_xy(0) = do2d_xy_time_count(0) + CEILING(             &
2232                                 ( end_time - MAX( skip_time_do2d_xy,         &
2233                                                   simulated_time_at_begin )  &
2234                                 ) / dt_do2d_xy )
2235             ELSE
2236                ntdim_2d_xy(1) = do2d_xy_time_count(1) + CEILING(             &
2237                                 ( end_time - MAX( skip_time_data_output_av,  &
2238                                                   simulated_time_at_begin )  &
2239                                 ) / dt_data_output_av )
2240             ENDIF
2241
2242!
2243!--          Check if the needed number of output time levels is increased
2244!--          compared to the number of time levels in the existing file.
[1308]2245             IF ( ntdim_2d_xy(av) > ntime_count )  THEN
2246                message_string = 'netCDF file for cross sections ' //          &
2247                                 TRIM( var ) // ' from previous run found,' // &
2248                                 '&but this file cannot be extended becaus' // &
2249                                 'e the number of output time levels&has b' // &
2250                                 'een increased compared to the previous s' // &
[1745]2251                                 'imulation.' //                               &
[1308]2252                                 '&New file is created instead.'
2253                CALL message( 'define_netcdf_header', 'PA0389', 0, 1, 0, 6, 0 )
2254                do2d_xy_time_count(av) = 0
2255                extend = .FALSE.
[1745]2256!
2257!--             Recalculate the needed time levels for the new file.
2258                IF ( av == 0 )  THEN
2259                   ntdim_2d_xy(0) = CEILING(                            &
2260                           ( end_time - MAX( skip_time_do2d_xy,         &
2261                                             simulated_time_at_begin )  &
2262                           ) / dt_do2d_xy )
2263                ELSE
2264                   ntdim_2d_xy(1) = CEILING(                            &
2265                           ( end_time - MAX( skip_time_data_output_av,  &
2266                                             simulated_time_at_begin )  &
2267                           ) / dt_data_output_av )
2268                ENDIF
[1308]2269                RETURN
2270             ENDIF
2271          ENDIF
2272
[1]2273!
2274!--       Dataset seems to be extendable.
2275!--       Now get the variable ids.
2276          i = 1
2277          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2278             IF ( INDEX( do2d(av,i), 'xy' ) /= 0 )  THEN
[519]2279                nc_stat = NF90_INQ_VARID( id_set_xy(av), do2d(av,i), &
[1]2280                                          id_var_do2d(av,i) )
[1783]2281                CALL netcdf_handle_error( 'netcdf_define_header', 138 )
[1031]2282#if defined( __netcdf4_parallel )
[493]2283!
2284!--             Set collective io operations for parallel io
[1031]2285                IF ( netcdf_data_format > 4 )  THEN
[493]2286                   nc_stat = NF90_VAR_PAR_ACCESS( id_set_xy(av),     &
2287                                                  id_var_do2d(av,i), &
2288                                                  NF90_COLLECTIVE )
[1783]2289                   CALL netcdf_handle_error( 'netcdf_define_header', 454 )
[493]2290                ENDIF
2291#endif
[1]2292             ENDIF
2293             i = i + 1
2294          ENDDO
2295
2296!
[359]2297!--       Update the title attribute on file
2298!--       In order to avoid 'data mode' errors if updated attributes are larger
2299!--       than their original size, NF90_PUT_ATT is called in 'define mode'
2300!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
2301!--       performance loss due to data copying; an alternative strategy would be
2302!--       to ensure equal attribute size in a job chain. Maybe revise later.
2303          IF ( av == 0 )  THEN
2304             time_average_text = ' '
2305          ELSE
2306             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
2307                                                            averaging_interval
2308          ENDIF
2309          nc_stat = NF90_REDEF( id_set_xy(av) )
[1783]2310          CALL netcdf_handle_error( 'netcdf_define_header', 431 )
[1]2311          nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'title', &
[359]2312                                  TRIM( run_description_header ) //    &
2313                                  TRIM( time_average_text ) )
[1783]2314          CALL netcdf_handle_error( 'netcdf_define_header', 139 )
[359]2315          nc_stat = NF90_ENDDEF( id_set_xy(av) )
[1783]2316          CALL netcdf_handle_error( 'netcdf_define_header', 432 )
[1031]2317          message_string = 'netCDF file for cross-sections ' //           &
[257]2318                            TRIM( var ) // ' from previous run found.' // &
2319                           '&This file will be extended.'
2320          CALL message( 'define_netcdf_header', 'PA0253', 0, 0, 0, 6, 0 )
2321         
[1]2322
2323       CASE ( 'xz_new' )
2324
2325!
2326!--       Define some global attributes of the dataset
2327          nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'Conventions', &
2328                                  'COARDS' )
[1783]2329          CALL netcdf_handle_error( 'netcdf_define_header', 140 )
[1]2330
2331          IF ( av == 0 )  THEN
2332             time_average_text = ' '
2333          ELSE
2334             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
2335                                                            averaging_interval
2336          ENDIF
2337          nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'title', &
2338                                  TRIM( run_description_header )  //   &
2339                                  TRIM( time_average_text ) )
[1783]2340          CALL netcdf_handle_error( 'netcdf_define_header', 141 )
[1]2341          IF ( av == 1 )  THEN
2342             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
2343             nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'time_avg', &
2344                                     TRIM( time_average_text ) )
[1783]2345             CALL netcdf_handle_error( 'netcdf_define_header', 141 )
[1]2346          ENDIF
2347
2348!
[1308]2349!--       Define time coordinate for xz sections.
2350!--       For parallel output the time dimensions has to be limited, otherwise
2351!--       the performance drops significantly.
2352          IF ( netcdf_data_format < 5 )  THEN
[1783]2353             CALL netcdf_create_dim( id_set_xz(av), 'time', NF90_UNLIMITED,    &
2354                                     id_dim_time_xz(av), 142 )
[1308]2355          ELSE
[1783]2356             CALL netcdf_create_dim( id_set_xz(av), 'time', ntdim_2d_xz(av),   &
2357                                     id_dim_time_xz(av), 525 )
[1308]2358          ENDIF
[1]2359
[1783]2360          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_time_xz(av) /),     &
2361                                  'time', NF90_DOUBLE, id_var_time_xz(av),     &
2362                                  'seconds', '', 143, 144, 000 )
[1]2363!
2364!--       Define the spatial dimensions and coordinates for xz-sections.
2365!--       First, determine the number of vertical sections to be written.
2366          IF ( section(1,2) == -9999 )  THEN
2367             RETURN
2368          ELSE
2369             ns = 1
2370             DO WHILE ( section(ns,2) /= -9999  .AND.  ns <= 100 )
2371                ns = ns + 1
2372             ENDDO
2373             ns = ns - 1
2374          ENDIF
2375
2376!
2377!--       Define y-axis (for scalar position)
[1783]2378          CALL netcdf_create_dim( id_set_xz(av), 'y_xz', ns, id_dim_y_xz(av),  &
2379                                  145 )
2380          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_y_xz(av) /),        &
2381                                  'y_xz', NF90_DOUBLE, id_var_y_xz(av),        &
2382                                  'meters', '', 146, 147, 000 )
[1]2383!
2384!--       Define y-axis (for v position)
[1783]2385          CALL netcdf_create_dim( id_set_xz(av), 'yv_xz', ns,                  &
2386                                  id_dim_yv_xz(av), 369 )
2387          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_yv_xz(av) /),       &
2388                                  'yv_xz', NF90_DOUBLE, id_var_yv_xz(av),      &
2389                                  'meters', '', 370, 371, 000 )
[1]2390!
2391!--       Define a variable to store the layer indices of the vertical cross
2392!--       sections
[1783]2393          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_y_xz(av) /),        &
2394                                  'ind_y_xz', NF90_DOUBLE,                     &
2395                                  id_var_ind_y_xz(av), 'gridpoints', '', 148,  &
2396                                  149, 000 )
[1]2397!
2398!--       Define x-axis (for scalar position)
[1783]2399          CALL netcdf_create_dim( id_set_xz(av), 'x', nx+2, id_dim_x_xz(av),   &
2400                                  150 )
2401          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_x_xz(av) /), 'x',   &
2402                                  NF90_DOUBLE, id_var_x_xz(av), 'meters', '',  &
2403                                  151, 152, 000 )
[1]2404!
2405!--       Define x-axis (for u position)
[1783]2406          CALL netcdf_create_dim( id_set_xz(av), 'xu', nx+2, id_dim_xu_xz(av), &
2407                                  372 )
2408          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_xu_xz(av) /), 'xu', &
2409                                  NF90_DOUBLE, id_var_xu_xz(av), 'meters', '', &
2410                                  373, 374, 000 )
[1]2411!
[1551]2412!--       Define the three z-axes (zu, zw, and zs)
[1783]2413          CALL netcdf_create_dim( id_set_xz(av), 'zu', nz+2, id_dim_zu_xz(av), &
2414                                  153 )
2415          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zu_xz(av) /), 'zu', &
2416                                  NF90_DOUBLE, id_var_zu_xz(av), 'meters', '', &
2417                                  154, 155, 000 )
2418          CALL netcdf_create_dim( id_set_xz(av), 'zw', nz+2, id_dim_zw_xz(av), &
2419                                  156 )
2420          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zw_xz(av) /), 'zw', &
2421                                  NF90_DOUBLE, id_var_zw_xz(av), 'meters', '', &
2422                                  157, 158, 000 )
[1]2423
[1551]2424          IF ( land_surface )  THEN
[1]2425
[1783]2426             CALL netcdf_create_dim( id_set_xz(av), 'zs', nzs,                 &
2427                                     id_dim_zs_xz(av), 542 )
2428             CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zs_xz(av) /),    &
2429                                     'zs', NF90_DOUBLE, id_var_zs_xz(av),      &
2430                                     'meters', '', 543, 544, 000 )
[1551]2431
2432          ENDIF
2433
[1]2434!
2435!--       Define the variables
2436          var_list = ';'
2437          i = 1
2438
2439          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2440
2441             IF ( INDEX( do2d(av,i), 'xz' ) /= 0 )  THEN
2442
2443!
2444!--             Check for the grid
2445                found = .TRUE.
2446                SELECT CASE ( do2d(av,i) )
2447!
2448!--                Most variables are defined on the zu grid
[1053]2449                   CASE ( 'e_xz', 'lpt_xz', 'nr_xz', 'p_xz', 'pc_xz', 'pr_xz', &
2450                          'prr_xz', 'pt_xz', 'q_xz', 'qc_xz', 'ql_xz',         &
2451                          'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz', 'qv_xz',  &
[1691]2452                          'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',                   &
2453                          'rad_sw_cs_hr_xz', 'rad_sw_hr_xz''rho_xz', 's_xz',   &
2454                          'sa_xz', 'vpt_xz' )
[1]2455
2456                      grid_x = 'x'
2457                      grid_y = 'y'
2458                      grid_z = 'zu'
2459!
2460!--                u grid
2461                   CASE ( 'u_xz' )
2462
2463                      grid_x = 'xu'
2464                      grid_y = 'y'
2465                      grid_z = 'zu'
2466!
2467!--                v grid
2468                   CASE ( 'v_xz' )
2469
2470                      grid_x = 'x'
2471                      grid_y = 'yv'
2472                      grid_z = 'zu'
2473!
2474!--                w grid
[1691]2475                   CASE ( 'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz',     &
2476                          'rad_sw_out_xz', 'w_xz' )
[1]2477
2478                      grid_x = 'x'
2479                      grid_y = 'y'
2480                      grid_z = 'zw'
2481
[1551]2482!
2483!--                soil grid
[1691]2484                   CASE ( 'm_soil_xz', 't_soil_xz' )
[1551]2485
2486                      grid_x = 'x'
2487                      grid_y = 'y'
2488                      grid_z = 'zs'
2489
[1]2490                   CASE DEFAULT
2491!
2492!--                   Check for user-defined quantities
2493                      CALL user_define_netcdf_grid( do2d(av,i), found, &
2494                                                    grid_x, grid_y, grid_z )
[1783]2495                      IF ( .NOT. found )  THEN
2496                         WRITE ( message_string, * ) 'no grid defined for', &
2497                                                ' variable ', TRIM( do2d(av,i) )
2498                         CALL message( 'define_netcdf_header', 'PA0244', &
2499                                       0, 1, 0, 6, 0 )
2500                      ENDIF
[1]2501
2502                END SELECT
2503
2504!
2505!--             Select the respective dimension ids
2506                IF ( grid_x == 'x' )  THEN
2507                   id_x = id_dim_x_xz(av)
2508                ELSEIF ( grid_x == 'xu' )  THEN
2509                   id_x = id_dim_xu_xz(av)
2510                ENDIF
2511
2512                IF ( grid_y == 'y' )  THEN
2513                   id_y = id_dim_y_xz(av)
2514                ELSEIF ( grid_y == 'yv' )  THEN
2515                   id_y = id_dim_yv_xz(av)
2516                ENDIF
2517
2518                IF ( grid_z == 'zu' )  THEN
2519                   id_z = id_dim_zu_xz(av)
2520                ELSEIF ( grid_z == 'zw' )  THEN
2521                   id_z = id_dim_zw_xz(av)
[1551]2522                ELSEIF ( grid_z == 'zs' )  THEN
2523                   id_z = id_dim_zs_xz(av)
[1]2524                ENDIF
2525
2526!
2527!--             Define the grid
[1783]2528                CALL netcdf_create_var( id_set_xz(av), (/ id_x, id_y, id_z,    &
2529                                        id_dim_time_xz(av) /), do2d(av,i),     &
2530                                        nc_precision(2), id_var_do2d(av,i),    &
2531                                        TRIM( do2d_unit(av,i) ), do2d(av,i),   &
2532                                        159, 160, 355 )
[1]2533
[1031]2534#if defined( __netcdf4_parallel )
[1308]2535
2536                IF ( netcdf_data_format > 4 )  THEN
[493]2537!
[1308]2538!--                Set no fill for every variable to increase performance.
2539                   nc_stat = NF90_DEF_VAR_FILL( id_set_xz(av),     &
2540                                                id_var_do2d(av,i), &
2541                                                1, 0 )
[1783]2542                   CALL netcdf_handle_error( 'netcdf_define_header', 534 )
[1308]2543!
2544!--                Set independent io operations for parallel io. Collective io
2545!--                is only allowed in case of a 1d-decomposition along x,
2546!--                because otherwise, not all PEs have output data.
[493]2547                   IF ( npey == 1 )  THEN
2548                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
2549                                                     id_var_do2d(av,i), &
2550                                                     NF90_COLLECTIVE )
2551                   ELSE
2552!
[1308]2553!--                   Test simulations showed that the output of cross sections
2554!--                   by all PEs in data_output_2d using NF90_COLLECTIVE is
2555!--                   faster than the output by the first row of PEs in
2556!--                   x-direction using NF90_INDEPENDENT.
2557                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),    & 
2558                                                    id_var_do2d(av,i), &
2559                                                    NF90_COLLECTIVE )
[493]2560!                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
2561!                                                     id_var_do2d(av,i), &
2562!                                                     NF90_INDEPENDENT )
2563                   ENDIF
[1783]2564                   CALL netcdf_handle_error( 'netcdf_define_header', 449 )
[493]2565                ENDIF
2566#endif
[1783]2567                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
2568
[1]2569             ENDIF
2570
2571             i = i + 1
2572
2573          ENDDO
2574
2575!
2576!--       No arrays to output. Close the netcdf file and return.
2577          IF ( i == 1 )  RETURN
2578
2579!
2580!--       Write the list of variables as global attribute (this is used by
2581!--       restart runs and by combine_plot_fields)
2582          nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'VAR_LIST', &
2583                                  var_list )
[1783]2584          CALL netcdf_handle_error( 'netcdf_define_header', 161 )
[1]2585
2586!
[1308]2587!--       Set general no fill, otherwise the performance drops significantly for
2588!--       parallel output.
2589          nc_stat = NF90_SET_FILL( id_set_xz(av), NF90_NOFILL, oldmode )
[1783]2590          CALL netcdf_handle_error( 'netcdf_define_header', 530 )
[1308]2591
2592!
[1031]2593!--       Leave netCDF define mode
[1]2594          nc_stat = NF90_ENDDEF( id_set_xz(av) )
[1783]2595          CALL netcdf_handle_error( 'netcdf_define_header', 162 )
[1]2596
2597!
[1308]2598!--       These data are only written by PE0 for parallel output to increase
2599!--       the performance.
2600          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
[1]2601
2602!
[1308]2603!--          Write axis data: y_xz, x, zu, zw
2604             ALLOCATE( netcdf_data(1:ns) )
[1]2605
2606!
[1308]2607!--          Write y_xz data (shifted by +dy/2)
2608             DO  i = 1, ns
2609                IF( section(i,2) == -1 )  THEN
[1353]2610                   netcdf_data(i) = -1.0_wp  ! section averaged along y
[1308]2611                ELSE
[1353]2612                   netcdf_data(i) = ( section(i,2) + 0.5_wp ) * dy
[1308]2613                ENDIF
2614             ENDDO
2615             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_y_xz(av), &
2616                                     netcdf_data, start = (/ 1 /),   &
2617                                     count = (/ ns /) )
[1783]2618             CALL netcdf_handle_error( 'netcdf_define_header', 163 )
[1]2619
2620!
[1308]2621!--          Write yv_xz data
2622             DO  i = 1, ns
2623                IF( section(i,2) == -1 )  THEN
[1353]2624                   netcdf_data(i) = -1.0_wp  ! section averaged along y
[1308]2625                ELSE
2626                   netcdf_data(i) = section(i,2) * dy
2627                ENDIF
2628             ENDDO
2629             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_yv_xz(av), &
2630                                     netcdf_data, start = (/ 1 /),    &
2631                                     count = (/ ns /) )
[1783]2632             CALL netcdf_handle_error( 'netcdf_define_header', 375 )
[1]2633
[1308]2634!
2635!--          Write gridpoint number data
2636             netcdf_data(1:ns) = section(1:ns,2)
2637             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_ind_y_xz(av), &
2638                                     netcdf_data, start = (/ 1 /),       &
2639                                     count = (/ ns /) )
[1783]2640             CALL netcdf_handle_error( 'netcdf_define_header', 164 )
[263]2641
[1]2642
[1308]2643             DEALLOCATE( netcdf_data )
2644
[1]2645!
[1308]2646!--          Write data for x (shifted by +dx/2) and xu axis
2647             ALLOCATE( netcdf_data(0:nx+1) )
[1]2648
[1308]2649             DO  i = 0, nx+1
[1353]2650                netcdf_data(i) = ( i + 0.5_wp ) * dx
[1308]2651             ENDDO
[1]2652
[1308]2653             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_x_xz(av), &
2654                                     netcdf_data, start = (/ 1 /),   &
2655                                     count = (/ nx+2 /) )
[1783]2656             CALL netcdf_handle_error( 'netcdf_define_header', 165 )
[1]2657
[1308]2658             DO  i = 0, nx+1
2659                netcdf_data(i) = i * dx
2660             ENDDO
[1]2661
[1308]2662             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_xu_xz(av), &
2663                                     netcdf_data, start = (/ 1 /),    &
2664                                     count = (/ nx+2 /) )
[1783]2665             CALL netcdf_handle_error( 'netcdf_define_header', 377 )
[1]2666
[1308]2667             DEALLOCATE( netcdf_data )
[1]2668
2669!
[1308]2670!--          Write zu and zw data (vertical axes)
2671             ALLOCATE( netcdf_data(0:nz+1) )
[1]2672
[1308]2673             netcdf_data(0:nz+1) = zu(nzb:nzt+1)
2674             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zu_xz(av), &
2675                                     netcdf_data, start = (/ 1 /),    &
2676                                     count = (/ nz+2 /) )
[1783]2677             CALL netcdf_handle_error( 'netcdf_define_header', 166 )
[1]2678
[1308]2679             netcdf_data(0:nz+1) = zw(nzb:nzt+1)
2680             nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zw_xz(av), &
2681                                     netcdf_data, start = (/ 1 /),    &
2682                                     count = (/ nz+2 /) )
[1783]2683             CALL netcdf_handle_error( 'netcdf_define_header', 167 )
[1]2684
[1551]2685!
2686!--          Write zs data
2687             IF ( land_surface )  THEN
2688                netcdf_data(0:nzs-1) = - zs(nzb_soil:nzt_soil)
2689                nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zs_xz(av), &
2690                                        netcdf_data(0:nzs), start = (/ 1 /),    &
2691                                        count = (/ nzt_soil-nzb_soil+1 /) )
[1783]2692               CALL netcdf_handle_error( 'netcdf_define_header', 548 )
[1551]2693             ENDIF
2694
2695
[1308]2696             DEALLOCATE( netcdf_data )
[1]2697
[1308]2698          ENDIF
[1]2699
[1308]2700
[1]2701       CASE ( 'xz_ext' )
2702
2703!
2704!--       Get the list of variables and compare with the actual run.
2705!--       First var_list_old has to be reset, since GET_ATT does not assign
2706!--       trailing blanks.
2707          var_list_old = ' '
2708          nc_stat = NF90_GET_ATT( id_set_xz(av), NF90_GLOBAL, 'VAR_LIST', &
2709                                  var_list_old )
[1783]2710          CALL netcdf_handle_error( 'netcdf_define_header', 168 )
[1]2711
2712          var_list = ';'
2713          i = 1
2714          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2715             IF ( INDEX( do2d(av,i), 'xz' ) /= 0 )  THEN
[519]2716                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
[1]2717             ENDIF
2718             i = i + 1
2719          ENDDO
2720
2721          IF ( av == 0 )  THEN
2722             var = '(xz)'
2723          ELSE
2724             var = '(xz_av)'
2725          ENDIF
2726
2727          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
[1031]2728             message_string = 'netCDF file for cross-sections ' //           &
[292]2729                              TRIM( var ) // ' from previous run found,' //  &
[257]2730                              '& but this file cannot be extended due to' // &
2731                              ' variable mismatch.' //                       &
2732                              '&New file is created instead.'
2733             CALL message( 'define_netcdf_header', 'PA0249', 0, 1, 0, 6, 0 )
[1]2734             extend = .FALSE.
2735             RETURN
2736          ENDIF
2737
2738!
2739!--       Calculate the number of current sections
2740          ns = 1
2741          DO WHILE ( section(ns,2) /= -9999  .AND.  ns <= 100 )
2742             ns = ns + 1
2743          ENDDO
2744          ns = ns - 1
2745
2746!
2747!--       Get and compare the number of vertical cross sections
2748          nc_stat = NF90_INQ_VARID( id_set_xz(av), 'y_xz', id_var_y_xz(av) )
[1783]2749          CALL netcdf_handle_error( 'netcdf_define_header', 169 )
[1]2750
2751          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_y_xz(av), &
2752                                           dimids = id_dim_y_xz_old )
[1783]2753          CALL netcdf_handle_error( 'netcdf_define_header', 170 )
[1]2754          id_dim_y_xz(av) = id_dim_y_xz_old(1)
2755
2756          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_y_xz(av), &
2757                                            len = ns_old )
[1783]2758          CALL netcdf_handle_error( 'netcdf_define_header', 171 )
[1]2759
2760          IF ( ns /= ns_old )  THEN
[1031]2761             message_string = 'netCDF file for cross-sections ' //          &
[292]2762                              TRIM( var ) // ' from previous run found,' // &
[257]2763                              '&but this file cannot be extended due to' // &
2764                              ' mismatch in number of' //                   & 
2765                              '&cross sections.' //                         &
2766                              '&New file is created instead.'
2767             CALL message( 'define_netcdf_header', 'PA0250', 0, 1, 0, 6, 0 )
[1]2768             extend = .FALSE.
2769             RETURN
2770          ENDIF
2771
2772!
2773!--       Get and compare the heights of the cross sections
2774          ALLOCATE( netcdf_data(1:ns_old) )
2775
2776          nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_y_xz(av), netcdf_data )
[1783]2777          CALL netcdf_handle_error( 'netcdf_define_header', 172 )
[1]2778
2779          DO  i = 1, ns
2780             IF ( section(i,2) /= -1 )  THEN
[600]2781                IF ( ( ( section(i,2) + 0.5 ) * dy ) /= netcdf_data(i) )  THEN
[1031]2782                   message_string = 'netCDF file for cross-sections ' //     &
[292]2783                               TRIM( var ) // ' from previous run found,' // &
[274]2784                               '&but this file cannot be extended' //        &
2785                               ' due to mismatch in cross' //                &
2786                               '&section levels.' //                         &
2787                                '&New file is created instead.'
2788                   CALL message( 'define_netcdf_header', 'PA0251', &
2789                                                                 0, 1, 0, 6, 0 )
[1]2790                   extend = .FALSE.
2791                   RETURN
2792                ENDIF
2793             ELSE
[1353]2794                IF ( -1.0_wp /= netcdf_data(i) )  THEN
[1031]2795                   message_string = 'netCDF file for cross-sections ' //     &
[292]2796                               TRIM( var ) // ' from previous run found,' // &
[274]2797                               '&but this file cannot be extended' //        &
2798                               ' due to mismatch in cross' //                &
2799                               '&section levels.' //                         &
2800                               '&New file is created instead.'
2801                   CALL message( 'define_netcdf_header', 'PA0251', &
2802                                                                 0, 1, 0, 6, 0 )
[1]2803                   extend = .FALSE.
2804                   RETURN
2805                ENDIF
2806             ENDIF
2807          ENDDO
2808
2809          DEALLOCATE( netcdf_data )
2810
2811!
2812!--       Get the id of the time coordinate (unlimited coordinate) and its
2813!--       last index on the file. The next time level is do2d..count+1.
2814!--       The current time must be larger than the last output time
2815!--       on the file.
2816          nc_stat = NF90_INQ_VARID( id_set_xz(av), 'time', id_var_time_xz(av) )
[1783]2817          CALL netcdf_handle_error( 'netcdf_define_header', 173 )
[1]2818
2819          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_time_xz(av), &
2820                                           dimids = id_dim_time_old )
[1783]2821          CALL netcdf_handle_error( 'netcdf_define_header', 174 )
[1]2822          id_dim_time_xz(av) = id_dim_time_old(1)
2823
2824          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_time_xz(av), &
[1308]2825                                            len = ntime_count )
[1783]2826          CALL netcdf_handle_error( 'netcdf_define_header', 175 )
[1]2827
[1308]2828!
2829!--       For non-parallel output use the last output time level of the netcdf
2830!--       file because the time dimension is unlimited. In case of parallel
2831!--       output the variable ntime_count could get the value of 9*10E36 because
2832!--       the time dimension is limited.
2833          IF ( netcdf_data_format < 5 ) do2d_xz_time_count(av) = ntime_count
2834
[1]2835          nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_time_xz(av),    &
2836                                  last_time_coordinate,                 &
2837                                  start = (/ do2d_xz_time_count(av) /), &
2838                                  count = (/ 1 /) )
[1783]2839          CALL netcdf_handle_error( 'netcdf_define_header', 176 )
[1]2840
2841          IF ( last_time_coordinate(1) >= simulated_time )  THEN
[1031]2842             message_string = 'netCDF file for cross sections ' //          &
[292]2843                              TRIM( var ) // ' from previous run found,' // &
[257]2844                              '&but this file cannot be extended becaus' // &
2845                              'e the current output time' //                &
2846                              '&is less or equal than the last output t' // &
2847                              'ime on this file.' //                        &
2848                              '&New file is created instead.'
2849             CALL message( 'define_netcdf_header', 'PA0252', 0, 1, 0, 6, 0 )
[1]2850             do2d_xz_time_count(av) = 0
2851             extend = .FALSE.
2852             RETURN
2853          ENDIF
[1745]2854
[1308]2855          IF ( netcdf_data_format > 4 )  THEN
[1745]2856!
2857!--          Set needed time levels (ntdim) to
2858!--          saved time levels + to be saved time levels.
2859             IF ( av == 0 )  THEN
2860                ntdim_2d_xz(0) = do2d_xz_time_count(0) + CEILING(             &
2861                                 ( end_time - MAX( skip_time_do2d_xz,         &
2862                                                   simulated_time_at_begin )  &
2863                                 ) / dt_do2d_xz )
2864             ELSE
2865                ntdim_2d_xz(1) = do2d_xz_time_count(1) + CEILING(             &
2866                                 ( end_time - MAX( skip_time_data_output_av,  &
2867                                                   simulated_time_at_begin )  &
2868                                 ) / dt_data_output_av )
2869             ENDIF
2870
2871!
2872!--          Check if the needed number of output time levels is increased
2873!--          compared to the number of time levels in the existing file.
[1308]2874             IF ( ntdim_2d_xz(av) > ntime_count )  THEN
2875                message_string = 'netCDF file for cross sections ' // &
2876                                 TRIM( var ) // ' from previous run found,' // &
2877                                 '&but this file cannot be extended becaus' // &
2878                                 'e the number of output time levels&has b' // &
2879                                 'een increased compared to the previous s' // &
[1745]2880                                 'imulation.' //                               &
[1308]2881                                 '&New file is created instead.'
2882                CALL message( 'define_netcdf_header', 'PA0390', 0, 1, 0, 6, 0 )
2883                do2d_xz_time_count(av) = 0
2884                extend = .FALSE.
[1745]2885!
2886!--             Recalculate the needed time levels for the new file.
2887                IF ( av == 0 )  THEN
2888                   ntdim_2d_xz(0) = CEILING(                            &
2889                           ( end_time - MAX( skip_time_do2d_xz,         &
2890                                             simulated_time_at_begin )  &
2891                           ) / dt_do2d_xz )
2892                ELSE
2893                   ntdim_2d_xz(1) = CEILING(                            &
2894                           ( end_time - MAX( skip_time_data_output_av,  &
2895                                             simulated_time_at_begin )  &
2896                           ) / dt_data_output_av )
2897                ENDIF
[1308]2898                RETURN
2899             ENDIF
2900          ENDIF
[1]2901
2902!
2903!--       Dataset seems to be extendable.
2904!--       Now get the variable ids.
2905          i = 1
2906          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2907             IF ( INDEX( do2d(av,i), 'xz' ) /= 0 )  THEN
[519]2908                nc_stat = NF90_INQ_VARID( id_set_xz(av), do2d(av,i), &
[1]2909                                          id_var_do2d(av,i) )
[1783]2910                CALL netcdf_handle_error( 'netcdf_define_header', 177 )
[1031]2911#if defined( __netcdf4_parallel )
[493]2912!
2913!--             Set independent io operations for parallel io. Collective io
2914!--             is only allowed in case of a 1d-decomposition along x, because
2915!--             otherwise, not all PEs have output data.
[1031]2916                IF ( netcdf_data_format > 4 )  THEN
[493]2917                   IF ( npey == 1 )  THEN
2918                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
2919                                                     id_var_do2d(av,i), &
2920                                                     NF90_COLLECTIVE )
2921                   ELSE
2922!
[1308]2923!--                   Test simulations showed that the output of cross sections
2924!--                   by all PEs in data_output_2d using NF90_COLLECTIVE is
2925!--                   faster than the output by the first row of PEs in
2926!--                   x-direction using NF90_INDEPENDENT.
[493]2927                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
2928                                                     id_var_do2d(av,i), &
2929                                                     NF90_COLLECTIVE )
2930!                      nc_stat = NF90_VAR_PAR_ACCESS( id_set_xz(av),     &
2931!                                                     id_var_do2d(av,i), &
2932!                                                     NF90_INDEPENDENT )
2933                   ENDIF
[1783]2934                   CALL netcdf_handle_error( 'netcdf_define_header', 455 )
[493]2935                ENDIF
2936#endif
[1]2937             ENDIF
2938             i = i + 1
2939          ENDDO
2940
2941!
[359]2942!--       Update the title attribute on file
2943!--       In order to avoid 'data mode' errors if updated attributes are larger
2944!--       than their original size, NF90_PUT_ATT is called in 'define mode'
2945!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
2946!--       performance loss due to data copying; an alternative strategy would be
2947!--       to ensure equal attribute size in a job chain. Maybe revise later.
2948          IF ( av == 0 )  THEN
2949             time_average_text = ' '
2950          ELSE
2951             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
2952                                                            averaging_interval
2953          ENDIF
2954          nc_stat = NF90_REDEF( id_set_xz(av) )
[1783]2955          CALL netcdf_handle_error( 'netcdf_define_header', 433 )
[1]2956          nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'title', &
[359]2957                                  TRIM( run_description_header ) //    &
2958                                  TRIM( time_average_text ) )
[1783]2959          CALL netcdf_handle_error( 'netcdf_define_header', 178 )
[359]2960          nc_stat = NF90_ENDDEF( id_set_xz(av) )
[1783]2961          CALL netcdf_handle_error( 'netcdf_define_header', 434 )
[1031]2962          message_string = 'netCDF file for cross-sections ' //           &
[257]2963                            TRIM( var ) // ' from previous run found.' // &
2964                           '&This file will be extended.'
2965          CALL message( 'define_netcdf_header', 'PA0253', 0, 0, 0, 6, 0 )
[1]2966
2967
2968       CASE ( 'yz_new' )
2969
2970!
2971!--       Define some global attributes of the dataset
2972          nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'Conventions', &
2973                                  'COARDS' )
[1783]2974          CALL netcdf_handle_error( 'netcdf_define_header', 179 )
[1]2975
2976          IF ( av == 0 )  THEN
2977             time_average_text = ' '
2978          ELSE
2979             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
2980                                                            averaging_interval
2981          ENDIF
2982          nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'title', &
2983                                  TRIM( run_description_header ) //    &
2984                                  TRIM( time_average_text ) )
[1783]2985          CALL netcdf_handle_error( 'netcdf_define_header', 180 )
[1]2986          IF ( av == 1 )  THEN
2987             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
2988             nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'time_avg', &
2989                                     TRIM( time_average_text ) )
[1783]2990             CALL netcdf_handle_error( 'netcdf_define_header', 180 )
[1]2991          ENDIF
2992
2993!
[1308]2994!--       Define time coordinate for yz sections.
2995!--       For parallel output the time dimensions has to be limited, otherwise
2996!--       the performance drops significantly.
2997          IF ( netcdf_data_format < 5 )  THEN
[1783]2998             CALL netcdf_create_dim( id_set_yz(av), 'time', NF90_UNLIMITED,    &
2999                                     id_dim_time_yz(av), 181 )
[1308]3000          ELSE
[1783]3001             CALL netcdf_create_dim( id_set_yz(av), 'time', ntdim_2d_yz(av),   &
3002                                     id_dim_time_yz(av), 526 )
[1308]3003          ENDIF
[1]3004
[1783]3005          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_time_yz(av) /),     &
3006                                  'time', NF90_DOUBLE, id_var_time_yz(av),     &
3007                                  'seconds', '', 182, 183, 000 )
[1]3008!
3009!--       Define the spatial dimensions and coordinates for yz-sections.
3010!--       First, determine the number of vertical sections to be written.
3011          IF ( section(1,3) == -9999 )  THEN
3012             RETURN
3013          ELSE
3014             ns = 1
3015             DO WHILE ( section(ns,3) /= -9999  .AND.  ns <= 100 )
3016                ns = ns + 1
3017             ENDDO
3018             ns = ns - 1
3019          ENDIF
3020
3021!
3022!--       Define x axis (for scalar position)
[1783]3023          CALL netcdf_create_dim( id_set_yz(av), 'x_yz', ns, id_dim_x_yz(av),  &
3024                                  184 )
3025          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_x_yz(av) /),        &
3026                                  'x_yz', NF90_DOUBLE, id_var_x_yz(av),        &
3027                                  'meters', '', 185, 186, 000 )
[1]3028!
3029!--       Define x axis (for u position)
[1783]3030          CALL netcdf_create_dim( id_set_yz(av), 'xu_yz', ns,                  &
3031                                  id_dim_xu_yz(av), 377 )
3032          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_xu_yz(av) /),       &
3033                                  'xu_yz', NF90_DOUBLE, id_var_xu_yz(av),      &
3034                                  'meters', '', 378, 379, 000 )
[1]3035!
3036!--       Define a variable to store the layer indices of the vertical cross
3037!--       sections
[1783]3038          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_x_yz(av) /),        &
3039                                  'ind_x_yz', NF90_DOUBLE,                     &
3040                                  id_var_ind_x_yz(av), 'gridpoints', '', 187,  &
3041                                  188, 000 )
[1]3042!
3043!--       Define y-axis (for scalar position)
[1783]3044          CALL netcdf_create_dim( id_set_yz(av), 'y', ny+2, id_dim_y_yz(av),   &
3045                                  189 )
3046          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_y_yz(av) /), 'y',   &
3047                                  NF90_DOUBLE, id_var_y_yz(av), 'meters', '',  &
3048                                  190, 191, 000 )
[1]3049!
3050!--       Define y-axis (for v position)
[1783]3051          CALL netcdf_create_dim( id_set_yz(av), 'yv', ny+2, id_dim_yv_yz(av), &
3052                                  380 )
3053          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_yv_yz(av) /), 'yv', &
3054                                  NF90_DOUBLE, id_var_yv_yz(av), 'meters', '', &
3055                                  381, 382, 000 )
[1]3056!
3057!--       Define the two z-axes (zu and zw)
[1783]3058          CALL netcdf_create_dim( id_set_yz(av), 'zu', nz+2, id_dim_zu_yz(av), &
3059                                  192 )
3060          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_zu_yz(av) /), 'zu', &
3061                                  NF90_DOUBLE, id_var_zu_yz(av), 'meters', '', &
3062                                  193, 194, 000 )
[1]3063
[1783]3064          CALL netcdf_create_dim( id_set_yz(av),