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

Last change on this file since 1031 was 1031, checked in by raasch, 12 years ago

netCDF4 without parallel file support implemented
additional define string netcdf4_parallel is required to switch on parallel file support

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