source: palm/trunk/SOURCE/netcdf.f90 @ 1746

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

last commit documented

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