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

Last change on this file since 48 was 48, checked in by raasch, 14 years ago

preliminary version, several changes to be explained later

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