source: palm/tags/release-3.3/SOURCE/netcdf.f90 @ 708

Last change on this file since 708 was 98, checked in by raasch, 17 years ago

updating comments and rc-file

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