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

Last change on this file since 275 was 274, checked in by heinze, 16 years ago

Indentation of the message calls corrected

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