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

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