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

Last change on this file since 1308 was 1308, checked in by fricke, 10 years ago

Adjustments for parallel NetCDF output for lccrayh/lccrayb (Cray XC30 systems)

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