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

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