source: palm/tags/release-3.1c/SOURCE/netcdf.f90 @ 112

Last change on this file since 112 was 4, checked in by raasch, 17 years ago

Id keyword set as property for all *.f90 files

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