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

Last change on this file since 170 was 144, checked in by letzel, 17 years ago

User-defined spectra.

Bugfix: extra '*' removed in user_statistics sample code (user_interface).

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