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

Last change on this file since 249 was 226, checked in by raasch, 16 years ago

preparations for the next release

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