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

Last change on this file since 181 was 175, checked in by steinfeld, 16 years ago

NCL scripts for spectra added. Bug fix concerning spectra in netcdf.f90

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