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

Last change on this file since 205 was 198, checked in by raasch, 16 years ago

file headers updated for the next release 3.5

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