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

Last change on this file since 564 was 564, checked in by helmke, 11 years ago

several changes for an unlimited output of mask data and message IDs changed

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