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

Last change on this file since 216 was 216, checked in by raasch, 13 years ago

reading mechanism for restart files completely revised

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