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

Last change on this file since 590 was 565, checked in by helmke, 14 years ago

last commit documented

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