source: palm/tags/release-3.8/SOURCE/netcdf.f90 @ 4343

Last change on this file since 4343 was 601, checked in by raasch, 11 years ago

last commit documented

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