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

Last change on this file since 290 was 288, checked in by heinze, 16 years ago

NetCDF unit attribute in timeseries output added, typographical error in modules fixed

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