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

Last change on this file since 993 was 993, checked in by hoffmann, 12 years ago

last commit documented

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