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

Last change on this file since 1981 was 1981, checked in by suehring, 5 years ago

last commit documented

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