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

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