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

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