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

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

branch revision comments from Marcus (rev 410) replaced by normal revision comments

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