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

Last change on this file since 1 was 1, checked in by raasch, 15 years ago

Initial repository layout and content

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