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

Last change on this file since 992 was 992, checked in by hoffmann, 9 years ago

removal of two informative messages in netcdf.f90

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