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

Last change on this file since 1354 was 1354, checked in by heinze, 8 years ago

last commit documented

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