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

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

last commit documented

  • Property svn:keywords set to Id
File size: 185.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!
10!
11! Former revisions:
12! -----------------
13! $Id: netcdf.f90 984 2012-08-21 14:21:34Z 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 ::  cross_profiles_log, 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          cross_profiles_log = .FALSE.
3203          message_string = ''
3204
3205          DO  i = 1, cross_profiles_count
3206             DO  j = 1, dopr_n
3207
3208                IF ( TRIM(cross_profiles_char(i)) == TRIM(data_output_pr(j)) ) &
3209                THEN
3210                   EXIT
3211                ENDIF
3212
3213                IF ( j == dopr_n )  THEN
3214                   cross_profiles_numb(i) = 999999
3215                   cross_profiles_log = .TRUE.
3216                   message_string = TRIM( message_string ) // ', ' // &
3217                                    TRIM( cross_profiles_char(i) )
3218                ENDIF
3219
3220             ENDDO
3221          ENDDO
3222
3223          IF ( cross_profiles_log )  THEN
3224             message_string = TRIM( message_string ) // ' is/are not' // &
3225                              ' defined in data_output_pr.'
3226             CALL message( 'define_netcdf_header', 'PA0352', 0, 0, 0, 6, 0 )
3227          ENDIF
3228
3229          DO i = 1, crmax
3230             IF ( cross_profiles_numb(i) == 999999 ) THEN
3231                DO j = i + 1, crmax
3232                   IF ( cross_profiles_numb(j) /= 999999 ) THEN
3233                      cross_profiles_char(i) = cross_profiles_char(j)
3234                      cross_profiles_numb(i) = cross_profiles_numb(j)
3235                      cross_profiles_numb(j) = 999999
3236                      EXIT
3237                   ENDIF
3238                ENDDO
3239             ENDIF
3240          ENDDO
3241
3242          DO i = 1, crmax-1
3243             IF ( cross_profiles_numb(i + 1) == 999999 ) THEN
3244                cross_profiles_count = i
3245                EXIT
3246             ENDIF
3247          ENDDO
3248!
3249!--       Check if all profiles defined in data_output_pr are defined in
3250!--       cross_profiles. If not, they will be added to cross_profiles.
3251          cross_profiles_log = .FALSE.
3252          message_string = ' '
3253
3254          DO  i = 1, dopr_n
3255             DO  j = 1, cross_profiles_count
3256
3257                IF ( TRIM(cross_profiles_char(j)) == TRIM(data_output_pr(i))) &
3258                THEN
3259                   EXIT
3260                ENDIF
3261
3262                IF (( j == cross_profiles_count ) .AND.                        &
3263                    ( cross_profiles_count <= crmax - 1))  THEN
3264                   cross_profiles_count = cross_profiles_count + 1
3265                   cross_profiles_maxi  = cross_profiles_maxi  + 1
3266                   cross_profiles_char(MIN( crmax, cross_profiles_count )) =  &
3267                                                      TRIM( data_output_pr(i) )
3268                   cross_profiles_numb(MIN( crmax, cross_profiles_count )) =  &
3269                                                      cross_profiles_maxi
3270                   cross_profiles_log   = .TRUE.
3271                   message_string = TRIM( message_string ) // ', ' // &
3272                                    TRIM( data_output_pr(i) )
3273                ENDIF
3274
3275             ENDDO
3276          ENDDO
3277
3278          IF ( cross_profiles_log )  THEN
3279             message_string = TRIM(message_string(2:)) // ' has/have been' //  &
3280                              ' added to cross_profiles.'
3281             CALL message( 'define_netcdf_header', 'PA0353', 0, 0, 0, 6, 0 )
3282          ENDIF
3283
3284          IF ( cross_profiles_count >= crmax )  THEN
3285             message_string = 'It is not allowed to arrange more than '&
3286                              // '100 profiles with cross_profiles. Apart '&
3287                              // 'from that, all profiles are saved to the '&
3288                              // 'netCDF file.'
3289             CALL message( 'define_netcdf_header', 'PA0354', 0, 0, 0, 6, 0 )
3290          ENDIF
3291
3292!
3293!--       Writing cross_profiles to netcdf header. This information can be
3294!--       used by palmplot. Each profile is separated by ",", each cross is
3295!--       separated by ";".
3296          char_cross_profiles = ';'
3297          id_last = 1
3298          cross_profiles_count = MIN( cross_profiles_count, crmax )
3299
3300          DO  i = 1, cross_profiles_count
3301
3302             IF ( cross_profiles_numb(i) /= 999999 )  THEN
3303                IF ( TRIM( char_cross_profiles ) == ';' )  THEN
3304                   char_cross_profiles = TRIM( char_cross_profiles ) // &
3305                                         TRIM( cross_profiles_char(i) )
3306                ELSEIF ( id_last == cross_profiles_numb(i) )  THEN
3307                   char_cross_profiles = TRIM( char_cross_profiles ) // &
3308                                         ',' // TRIM( cross_profiles_char(i) )
3309                ELSE
3310                   char_cross_profiles = TRIM( char_cross_profiles ) // &
3311                                         ';' // TRIM( cross_profiles_char(i) )
3312                ENDIF
3313                id_last = cross_profiles_numb(i)
3314             ENDIF
3315
3316          ENDDO
3317
3318          char_cross_profiles = TRIM( char_cross_profiles ) // ';'
3319
3320          nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'cross_profiles',   &
3321                                  TRIM( char_cross_profiles ) )
3322          CALL handle_netcdf_error( 'netcdf', 521 )
3323
3324!
3325!--       Define time coordinate for profiles (unlimited dimension)
3326          nc_stat = NF90_DEF_DIM( id_set_pr, 'time', NF90_UNLIMITED, &
3327                                  id_dim_time_pr )
3328          CALL handle_netcdf_error( 'netcdf', 220 )
3329
3330          nc_stat = NF90_DEF_VAR( id_set_pr, 'time', NF90_DOUBLE, &
3331                                  id_dim_time_pr, id_var_time_pr )
3332          CALL handle_netcdf_error( 'netcdf', 221 )
3333
3334          nc_stat = NF90_PUT_ATT( id_set_pr, id_var_time_pr, 'units', 'seconds')
3335          CALL handle_netcdf_error( 'netcdf', 222 )
3336
3337!
3338!--       Define the variables
3339          var_list = ';'
3340          DO  i = 1, dopr_n
3341
3342             IF ( statistic_regions == 0 )  THEN
3343
3344!
3345!--             Define the z-axes (each variable gets its own z-axis)
3346                nc_stat = NF90_DEF_DIM( id_set_pr,                        &
3347                                        'z' // TRIM( data_output_pr(i) ), &
3348                                        nzt+2-nzb, id_dim_z_pr(i,0) )
3349                CALL handle_netcdf_error( 'netcdf', 223 )
3350
3351                nc_stat = NF90_DEF_VAR( id_set_pr,                        &
3352                                        'z' // TRIM( data_output_pr(i) ), &
3353                                        NF90_DOUBLE, id_dim_z_pr(i,0),    &
3354                                        id_var_z_pr(i,0) )
3355                CALL handle_netcdf_error( 'netcdf', 224 )
3356
3357                nc_stat = NF90_PUT_ATT( id_set_pr, id_var_z_pr(i,0), 'units', &
3358                                        'meters' )
3359                CALL handle_netcdf_error( 'netcdf', 225 )
3360!
3361!--             Define the variable
3362                nc_stat = NF90_DEF_VAR( id_set_pr, data_output_pr(i),         &
3363                                        nc_precision(5), (/ id_dim_z_pr(i,0), &
3364                                        id_dim_time_pr /), id_var_dopr(i,0) )
3365                CALL handle_netcdf_error( 'netcdf', 226 )
3366
3367                nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,0), &
3368                                        'long_name', TRIM( data_output_pr(i) ) )
3369                CALL handle_netcdf_error( 'netcdf', 227 )
3370
3371                nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,0), &
3372                                        'units', TRIM( dopr_unit(i) ) )
3373                CALL handle_netcdf_error( 'netcdf', 228 )
3374
3375                var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) //  ';'
3376
3377             ELSE
3378!
3379!--             If statistic regions are defined, add suffix _SR+#SR to the
3380!--             names
3381                DO  j = 0, statistic_regions
3382                   WRITE ( suffix, '(''_'',I1)' )  j
3383
3384!
3385!--                Define the z-axes (each variable gets it own z-axis)
3386                   nc_stat = NF90_DEF_DIM( id_set_pr,                          &
3387                                         'z'//TRIM(data_output_pr(i))//suffix, &
3388                                           nzt+2-nzb, id_dim_z_pr(i,j) )
3389                   CALL handle_netcdf_error( 'netcdf', 229 )
3390
3391                   nc_stat = NF90_DEF_VAR( id_set_pr,                          &
3392                                         'z'//TRIM(data_output_pr(i))//suffix, &
3393                                           nc_precision(5), id_dim_z_pr(i,j),  &
3394                                           id_var_z_pr(i,j) )
3395                   CALL handle_netcdf_error( 'netcdf', 230 )
3396
3397                   nc_stat = NF90_PUT_ATT( id_set_pr, id_var_z_pr(i,j), &
3398                                           'units', 'meters' )
3399                   CALL handle_netcdf_error( 'netcdf', 231 )
3400
3401!
3402!--                Define the variable
3403                   nc_stat = NF90_DEF_VAR( id_set_pr,                         &
3404                                           TRIM(data_output_pr(i)) // suffix, &
3405                                           nc_precision(5),                   &
3406                                           (/ id_dim_z_pr(i,j),               &
3407                                           id_dim_time_pr /), id_var_dopr(i,j) )
3408                   CALL handle_netcdf_error( 'netcdf', 232 )
3409
3410                   nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,j),        &
3411                                           'long_name',                        &
3412                                           TRIM( data_output_pr(i) ) // ' SR ' &
3413                                           // suffix(2:2) )
3414                   CALL handle_netcdf_error( 'netcdf', 233 )
3415
3416                   nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,j), &
3417                                           'units', TRIM( dopr_unit(i) ) )
3418                   CALL handle_netcdf_error( 'netcdf', 234 )
3419
3420                   var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) // &
3421                              suffix // ';'
3422
3423                ENDDO
3424
3425             ENDIF
3426
3427          ENDDO
3428
3429!
3430!--       Write the list of variables as global attribute (this is used by
3431!--       restart runs)
3432          nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'VAR_LIST', var_list )
3433          CALL handle_netcdf_error( 'netcdf', 235 )
3434
3435!
3436!--       Define normalization variables (as time series)
3437          DO  i = 1, dopr_norm_num
3438
3439             nc_stat = NF90_DEF_VAR( id_set_pr, 'NORM_' // &
3440                                     TRIM( dopr_norm_names(i) ), &
3441                                     nc_precision(5), (/ id_dim_time_pr /), &
3442                                     id_var_norm_dopr(i) )
3443             CALL handle_netcdf_error( 'netcdf', 236 )
3444
3445             nc_stat = NF90_PUT_ATT( id_set_pr, id_var_norm_dopr(i), &
3446                                     'long_name',                    &
3447                                     TRIM( dopr_norm_longnames(i) ) )
3448             CALL handle_netcdf_error( 'netcdf', 237 )
3449
3450          ENDDO
3451
3452!
3453!--       Leave NetCDF define mode
3454          nc_stat = NF90_ENDDEF( id_set_pr )
3455          CALL handle_netcdf_error( 'netcdf', 238 )
3456
3457!
3458!--       Write z-axes data
3459          DO  i = 1, dopr_n
3460             DO  j = 0, statistic_regions
3461
3462                nc_stat = NF90_PUT_VAR( id_set_pr, id_var_z_pr(i,j),      &
3463                                        hom(nzb:nzt+1,2,dopr_index(i),0), &
3464                                        start = (/ 1 /),                  &
3465                                        count = (/ nzt-nzb+2 /) )
3466                CALL handle_netcdf_error( 'netcdf', 239 )
3467
3468             ENDDO
3469          ENDDO
3470
3471
3472       CASE ( 'pr_ext' )
3473
3474!
3475!--       Get the list of variables and compare with the actual run.
3476!--       First var_list_old has to be reset, since GET_ATT does not assign
3477!--       trailing blanks.
3478          var_list_old = ' '
3479          nc_stat = NF90_GET_ATT( id_set_pr, NF90_GLOBAL, 'VAR_LIST', &
3480                                  var_list_old )
3481          CALL handle_netcdf_error( 'netcdf', 240 )
3482
3483          var_list = ';'
3484          DO  i = 1, dopr_n
3485
3486             IF ( statistic_regions == 0 )  THEN
3487                var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) // ';'
3488             ELSE
3489                DO  j = 0, statistic_regions
3490                   WRITE ( suffix, '(''_'',I1)' )  j
3491                   var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) // &
3492                              suffix // ';'
3493                ENDDO
3494             ENDIF
3495
3496          ENDDO
3497
3498          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
3499             message_string = 'NetCDF file for vertical profiles ' //        &
3500                              'from previous run found,' //                  &
3501                              '& but this file cannot be extended due to' // &
3502                              ' variable mismatch.' //                       &
3503                              '&New file is created instead.'
3504             CALL message( 'define_netcdf_header', 'PA0254', 0, 1, 0, 6, 0 )
3505             extend = .FALSE.
3506             RETURN
3507          ENDIF
3508
3509!
3510!--       Get the id of the time coordinate (unlimited coordinate) and its
3511!--       last index on the file. The next time level is dopr..count+1.
3512!--       The current time must be larger than the last output time
3513!--       on the file.
3514          nc_stat = NF90_INQ_VARID( id_set_pr, 'time', id_var_time_pr )
3515          CALL handle_netcdf_error( 'netcdf', 241 )
3516
3517          nc_stat = NF90_INQUIRE_VARIABLE( id_set_pr, id_var_time_pr, &
3518                                           dimids = id_dim_time_old )
3519          CALL handle_netcdf_error( 'netcdf', 242 )
3520          id_dim_time_pr = id_dim_time_old(1)
3521
3522          nc_stat = NF90_INQUIRE_DIMENSION( id_set_pr, id_dim_time_pr, &
3523                                            len = dopr_time_count )
3524          CALL handle_netcdf_error( 'netcdf', 243 )
3525
3526          nc_stat = NF90_GET_VAR( id_set_pr, id_var_time_pr,        &
3527                                  last_time_coordinate,             &
3528                                  start = (/ dopr_time_count /), &
3529                                  count = (/ 1 /) )
3530          CALL handle_netcdf_error( 'netcdf', 244 )
3531
3532          IF ( last_time_coordinate(1) >= simulated_time )  THEN
3533             message_string = 'NetCDF file for vertical profiles ' //       &
3534                              'from previous run found,' //                 &
3535                              '&but this file cannot be extended becaus' // &
3536                              'e the current output time' //                &
3537                              '&is less or equal than the last output t' // &
3538                              'ime on this file.' //                        &
3539                              '&New file is created instead.'
3540             CALL message( 'define_netcdf_header', 'PA0255', 0, 1, 0, 6, 0 )
3541             dopr_time_count = 0
3542             extend = .FALSE.
3543             RETURN
3544          ENDIF
3545
3546!
3547!--       Dataset seems to be extendable.
3548!--       Now get the variable ids.
3549          i = 1
3550          DO  i = 1, dopr_n
3551 
3552             IF ( statistic_regions == 0 )  THEN
3553                nc_stat = NF90_INQ_VARID( id_set_pr, data_output_pr(i), &
3554                                          id_var_dopr(i,0) )
3555                CALL handle_netcdf_error( 'netcdf', 245 )
3556             ELSE
3557                DO  j = 0, statistic_regions
3558                   WRITE ( suffix, '(''_'',I1)' )  j
3559                   netcdf_var_name = TRIM( data_output_pr(i) ) // suffix
3560                   nc_stat = NF90_INQ_VARID( id_set_pr, netcdf_var_name, &
3561                                             id_var_dopr(i,j) )
3562                   CALL handle_netcdf_error( 'netcdf', 246 )
3563                ENDDO
3564             ENDIF
3565
3566          ENDDO
3567
3568!
3569!--       Get ids of the normalization variables
3570          DO  i = 1, dopr_norm_num
3571             nc_stat = NF90_INQ_VARID( id_set_pr,                             &
3572                                       'NORM_' // TRIM( dopr_norm_names(i) ), &
3573                                       id_var_norm_dopr(i) )
3574             CALL handle_netcdf_error( 'netcdf', 247 )
3575          ENDDO
3576
3577!
3578!--       Update the title attribute on file
3579!--       In order to avoid 'data mode' errors if updated attributes are larger
3580!--       than their original size, NF90_PUT_ATT is called in 'define mode'
3581!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
3582!--       performance loss due to data copying; an alternative strategy would be
3583!--       to ensure equal attribute size in a job chain. Maybe revise later.
3584          IF ( averaging_interval_pr == 0.0 )  THEN
3585             time_average_text = ' '
3586          ELSE
3587             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
3588                                                            averaging_interval_pr
3589          ENDIF
3590          nc_stat = NF90_REDEF( id_set_pr )
3591          CALL handle_netcdf_error( 'netcdf', 437 )
3592          nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'title', &
3593                                  TRIM( run_description_header ) //    &
3594                                  TRIM( time_average_text ) )
3595          CALL handle_netcdf_error( 'netcdf', 248 )
3596
3597          nc_stat = NF90_ENDDEF( id_set_pr )
3598          CALL handle_netcdf_error( 'netcdf', 438 )
3599          message_string = 'NetCDF file for vertical profiles ' // & 
3600                           'from previous run found.' //           &
3601                           '&This file will be extended.'
3602          CALL message( 'define_netcdf_header', 'PA0256', 0, 0, 0, 6, 0 )
3603
3604
3605       CASE ( 'ts_new' )
3606
3607!
3608!--       Define some global attributes of the dataset
3609          nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'title', &
3610                                  TRIM( run_description_header ) )
3611          CALL handle_netcdf_error( 'netcdf', 249 )
3612
3613!
3614!--       Define time coordinate for time series (unlimited dimension)
3615          nc_stat = NF90_DEF_DIM( id_set_ts, 'time', NF90_UNLIMITED, &
3616                                  id_dim_time_ts )
3617          CALL handle_netcdf_error( 'netcdf', 250 )
3618
3619          nc_stat = NF90_DEF_VAR( id_set_ts, 'time', NF90_DOUBLE, &
3620                                  id_dim_time_ts, id_var_time_ts )
3621          CALL handle_netcdf_error( 'netcdf', 251 )
3622
3623          nc_stat = NF90_PUT_ATT( id_set_ts, id_var_time_ts, 'units', 'seconds')
3624          CALL handle_netcdf_error( 'netcdf', 252 )
3625
3626!
3627!--       Define the variables
3628          var_list = ';'
3629          DO  i = 1, dots_num
3630
3631             IF ( statistic_regions == 0 )  THEN
3632
3633                nc_stat = NF90_DEF_VAR( id_set_ts, dots_label(i),              &
3634                                        nc_precision(6), (/ id_dim_time_ts /), &
3635                                        id_var_dots(i,0) )
3636                CALL handle_netcdf_error( 'netcdf', 253 )
3637
3638                nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,0), &
3639                                        'long_name', TRIM( dots_label(i) ) )
3640                CALL handle_netcdf_error( 'netcdf', 254 )
3641
3642                nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,0), &
3643                                        'units', TRIM( dots_unit(i) ) )
3644                CALL handle_netcdf_error( 'netcdf', 255 )
3645
3646                var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // ';'
3647
3648             ELSE
3649!
3650!--             If statistic regions are defined, add suffix _SR+#SR to the
3651!--             names
3652                DO  j = 0, statistic_regions
3653                   WRITE ( suffix, '(''_'',I1)' )  j
3654
3655                   nc_stat = NF90_DEF_VAR( id_set_ts,                         &
3656                                           TRIM( dots_label(i) ) // suffix,   &
3657                                           nc_precision(6),                   &
3658                                           (/ id_dim_time_ts /),              &
3659                                           id_var_dots(i,j) )
3660                   CALL handle_netcdf_error( 'netcdf', 256 )
3661
3662                   nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,j),       &
3663                                           'long_name',                       &
3664                                           TRIM( dots_label(i) ) // ' SR ' // &
3665                                           suffix(2:2) )
3666                   CALL handle_netcdf_error( 'netcdf', 257 )
3667
3668                   nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,j), &
3669                                           'units', TRIM( dots_unit(i) ) )
3670                   CALL handle_netcdf_error( 'netcdf', 347 )
3671
3672                   var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // &
3673                              suffix // ';'
3674
3675                ENDDO
3676
3677             ENDIF
3678
3679          ENDDO
3680
3681!
3682!--       Write the list of variables as global attribute (this is used by
3683!--       restart runs)
3684          nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'VAR_LIST', var_list )
3685          CALL handle_netcdf_error( 'netcdf', 258 )
3686
3687!
3688!--       Leave NetCDF define mode
3689          nc_stat = NF90_ENDDEF( id_set_ts )
3690          CALL handle_netcdf_error( 'netcdf', 259 )
3691
3692
3693       CASE ( 'ts_ext' )
3694
3695!
3696!--       Get the list of variables and compare with the actual run.
3697!--       First var_list_old has to be reset, since GET_ATT does not assign
3698!--       trailing blanks.
3699          var_list_old = ' '
3700          nc_stat = NF90_GET_ATT( id_set_ts, NF90_GLOBAL, 'VAR_LIST', &
3701                                  var_list_old )
3702          CALL handle_netcdf_error( 'netcdf', 260 )
3703
3704          var_list = ';'
3705          i = 1
3706          DO  i = 1, dots_num
3707
3708             IF ( statistic_regions == 0 )  THEN
3709                var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // ';'
3710             ELSE
3711                DO  j = 0, statistic_regions
3712                   WRITE ( suffix, '(''_'',I1)' )  j
3713                   var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // &
3714                              suffix // ';'
3715                ENDDO
3716             ENDIF
3717
3718          ENDDO
3719
3720          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
3721             message_string = 'NetCDF file for time series ' //              &
3722                              'from previous run found,' //                  &
3723                              '& but this file cannot be extended due to' // &
3724                              ' variable mismatch.' //                       &
3725                              '&New file is created instead.'
3726             CALL message( 'define_netcdf_header', 'PA0257', 0, 1, 0, 6, 0 )
3727             extend = .FALSE.
3728             RETURN
3729          ENDIF
3730
3731!
3732!--       Get the id of the time coordinate (unlimited coordinate) and its
3733!--       last index on the file. The next time level is dots..count+1.
3734!--       The current time must be larger than the last output time
3735!--       on the file.
3736          nc_stat = NF90_INQ_VARID( id_set_ts, 'time', id_var_time_ts )
3737          CALL handle_netcdf_error( 'netcdf', 261 )
3738
3739          nc_stat = NF90_INQUIRE_VARIABLE( id_set_ts, id_var_time_ts, &
3740                                           dimids = id_dim_time_old )
3741          CALL handle_netcdf_error( 'netcdf', 262 )
3742          id_dim_time_ts = id_dim_time_old(1)
3743
3744          nc_stat = NF90_INQUIRE_DIMENSION( id_set_ts, id_dim_time_ts, &
3745                                            len = dots_time_count )
3746          CALL handle_netcdf_error( 'netcdf', 263 )
3747
3748          nc_stat = NF90_GET_VAR( id_set_ts, id_var_time_ts,        &
3749                                  last_time_coordinate,             &
3750                                  start = (/ dots_time_count /), &
3751                                  count = (/ 1 /) )
3752          CALL handle_netcdf_error( 'netcdf', 264 )
3753
3754          IF ( last_time_coordinate(1) >= simulated_time )  THEN
3755             message_string = 'NetCDF file for time series ' //             &
3756                              'from previous run found,' //                 &
3757                              '&but this file cannot be extended becaus' // &
3758                              'e the current output time' //                &
3759                              '&is less or equal than the last output t' // &
3760                              'ime on this file.' //                        &
3761                              '&New file is created instead.'
3762             CALL message( 'define_netcdf_header', 'PA0258', 0, 1, 0, 6, 0 )
3763             dots_time_count = 0
3764             extend = .FALSE.
3765             RETURN
3766          ENDIF
3767
3768!
3769!--       Dataset seems to be extendable.
3770!--       Now get the variable ids
3771          i = 1
3772          DO  i = 1, dots_num
3773 
3774             IF ( statistic_regions == 0 )  THEN
3775                nc_stat = NF90_INQ_VARID( id_set_ts, dots_label(i), &
3776                                          id_var_dots(i,0) )
3777                CALL handle_netcdf_error( 'netcdf', 265 )
3778             ELSE
3779                DO  j = 0, statistic_regions
3780                   WRITE ( suffix, '(''_'',I1)' )  j
3781                   netcdf_var_name = TRIM( dots_label(i) ) // suffix
3782                   nc_stat = NF90_INQ_VARID( id_set_ts, netcdf_var_name, &
3783                                             id_var_dots(i,j) )
3784                   CALL handle_netcdf_error( 'netcdf', 266 )
3785                ENDDO
3786             ENDIF
3787
3788          ENDDO
3789
3790!
3791!--       Update the title attribute on file
3792!--       In order to avoid 'data mode' errors if updated attributes are larger
3793!--       than their original size, NF90_PUT_ATT is called in 'define mode'
3794!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
3795!--       performance loss due to data copying; an alternative strategy would be
3796!--       to ensure equal attribute size in a job chain. Maybe revise later.
3797          nc_stat = NF90_REDEF( id_set_ts )
3798          CALL handle_netcdf_error( 'netcdf', 439 )
3799          nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'title', &
3800                                  TRIM( run_description_header ) )
3801          CALL handle_netcdf_error( 'netcdf', 267 )
3802          nc_stat = NF90_ENDDEF( id_set_ts )
3803          CALL handle_netcdf_error( 'netcdf', 440 )
3804          message_string = 'NetCDF file for time series ' // & 
3805                           'from previous run found.' //     &
3806                           '&This file will be extended.'
3807          CALL message( 'define_netcdf_header', 'PA0259', 0, 0, 0, 6, 0 )
3808
3809
3810       CASE ( 'sp_new' )
3811
3812!
3813!--       Define some global attributes of the dataset
3814          IF ( averaging_interval_sp /= 0.0 )  THEN
3815             WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
3816                                                            averaging_interval_sp
3817             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title',  &
3818                                     TRIM( run_description_header ) // &
3819                                     TRIM( time_average_text ) )
3820             CALL handle_netcdf_error( 'netcdf', 268 )
3821
3822             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval_sp
3823             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'time_avg', &
3824                                     TRIM( time_average_text ) )
3825          ELSE
3826             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
3827                                     TRIM( run_description_header ) )
3828          ENDIF
3829          CALL handle_netcdf_error( 'netcdf', 269 )
3830
3831!
3832!--       Define time coordinate for spectra (unlimited dimension)
3833          nc_stat = NF90_DEF_DIM( id_set_sp, 'time', NF90_UNLIMITED, &
3834                                  id_dim_time_sp )
3835          CALL handle_netcdf_error( 'netcdf', 270 )
3836
3837          nc_stat = NF90_DEF_VAR( id_set_sp, 'time', NF90_DOUBLE, &
3838                                  id_dim_time_sp, id_var_time_sp )
3839          CALL handle_netcdf_error( 'netcdf', 271 )
3840
3841          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_time_sp, 'units', 'seconds')
3842          CALL handle_netcdf_error( 'netcdf', 272 )
3843
3844!
3845!--       Define the spatial dimensions and coordinates for spectra.
3846!--       First, determine the number of vertical levels for which spectra
3847!--       are to be output.
3848          ns = 1
3849          DO WHILE ( comp_spectra_level(ns) /= 999999  .AND.  ns <= 100 )
3850             ns = ns + 1
3851          ENDDO
3852          ns = ns - 1
3853
3854!
3855!--       Define vertical coordinate grid (zu grid)
3856          nc_stat = NF90_DEF_DIM( id_set_sp, 'zu_sp', ns, id_dim_zu_sp )
3857          CALL handle_netcdf_error( 'netcdf', 273 )
3858
3859          nc_stat = NF90_DEF_VAR( id_set_sp, 'zu_sp', NF90_DOUBLE, &
3860                                  id_dim_zu_sp, id_var_zu_sp )
3861          CALL handle_netcdf_error( 'netcdf', 274 )
3862
3863          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_zu_sp, 'units', 'meters' )
3864          CALL handle_netcdf_error( 'netcdf', 275 )
3865
3866!
3867!--       Define vertical coordinate grid (zw grid)
3868          nc_stat = NF90_DEF_DIM( id_set_sp, 'zw_sp', ns, id_dim_zw_sp )
3869          CALL handle_netcdf_error( 'netcdf', 276 )
3870
3871          nc_stat = NF90_DEF_VAR( id_set_sp, 'zw_sp', NF90_DOUBLE, &
3872                                  id_dim_zw_sp, id_var_zw_sp )
3873          CALL handle_netcdf_error( 'netcdf', 277 )
3874
3875          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_zw_sp, 'units', 'meters' )
3876          CALL handle_netcdf_error( 'netcdf', 278 )
3877
3878!
3879!--       Define x-axis
3880          nc_stat = NF90_DEF_DIM( id_set_sp, 'k_x', nx/2, id_dim_x_sp )
3881          CALL handle_netcdf_error( 'netcdf', 279 )
3882
3883          nc_stat = NF90_DEF_VAR( id_set_sp, 'k_x', NF90_DOUBLE, id_dim_x_sp, &
3884                                  id_var_x_sp )
3885          CALL handle_netcdf_error( 'netcdf', 280 )
3886
3887          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_x_sp, 'units', 'm-1' )
3888          CALL handle_netcdf_error( 'netcdf', 281 )
3889
3890!
3891!--       Define y-axis
3892          nc_stat = NF90_DEF_DIM( id_set_sp, 'k_y', ny/2, id_dim_y_sp )
3893          CALL handle_netcdf_error( 'netcdf', 282 )
3894
3895          nc_stat = NF90_DEF_VAR( id_set_sp, 'k_y', NF90_DOUBLE, id_dim_y_sp, &
3896                                  id_var_y_sp )
3897          CALL handle_netcdf_error( 'netcdf', 283 )
3898
3899          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_y_sp, 'units', 'm-1' )
3900          CALL handle_netcdf_error( 'netcdf', 284 )
3901
3902!
3903!--       Define the variables
3904          var_list = ';'
3905          i = 1
3906          DO WHILE ( data_output_sp(i) /= ' '  .AND.  i <= 10 )
3907!
3908!--          First check for the vertical grid
3909             found = .TRUE.
3910             SELECT CASE ( data_output_sp(i) )
3911!
3912!--             Most variables are defined on the zu levels
3913                CASE ( 'e', 'lpt', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c',  &
3914                       'ql_v', 'ql_vp', 'qv', 'rho', 's', 'sa', 'u', 'v', 'vpt' )
3915
3916                   grid_z = 'zu'
3917!
3918!--             zw levels
3919                CASE ( 'w' )
3920
3921                   grid_z = 'zw'
3922
3923                CASE DEFAULT
3924!
3925!--                Check for user-defined quantities (found, grid_x and grid_y
3926!--                are dummies)
3927                   CALL user_define_netcdf_grid( data_output_sp(i), found,  &
3928                                                 grid_x, grid_y, grid_z )
3929
3930             END SELECT
3931
3932             IF ( INDEX( spectra_direction(i), 'x' ) /= 0 )  THEN
3933
3934!
3935!--             Define the variable
3936                netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
3937                IF ( TRIM( grid_z ) == 'zw' )  THEN
3938                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
3939                                           nc_precision(7), (/ id_dim_x_sp, &
3940                                           id_dim_zw_sp, id_dim_time_sp /), &
3941                                           id_var_dospx(i) )
3942                ELSE
3943                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
3944                                           nc_precision(7), (/ id_dim_x_sp, &
3945                                           id_dim_zu_sp, id_dim_time_sp /), &
3946                                           id_var_dospx(i) )
3947                ENDIF
3948                CALL handle_netcdf_error( 'netcdf', 285 )
3949
3950                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospx(i), &
3951                                        'long_name', netcdf_var_name )
3952                CALL handle_netcdf_error( 'netcdf', 286 )
3953
3954                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospx(i), &
3955                                        'units', 'unknown' )
3956                CALL handle_netcdf_error( 'netcdf', 287 )
3957
3958                var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
3959
3960             ENDIF
3961
3962             IF ( INDEX( spectra_direction(i), 'y' ) /= 0 )  THEN
3963
3964!
3965!--             Define the variable
3966                netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
3967                IF ( TRIM( grid_z ) == 'zw' )  THEN
3968                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
3969                                           nc_precision(7), (/ id_dim_y_sp, &
3970                                           id_dim_zw_sp, id_dim_time_sp /), &
3971                                           id_var_dospy(i) )
3972                ELSE
3973                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
3974                                           nc_precision(7), (/ id_dim_y_sp, &
3975                                           id_dim_zu_sp, id_dim_time_sp /), &
3976                                           id_var_dospy(i) )
3977                ENDIF
3978                CALL handle_netcdf_error( 'netcdf', 288 )
3979
3980                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospy(i), &
3981                                        'long_name', netcdf_var_name )
3982                CALL handle_netcdf_error( 'netcdf', 289 )
3983
3984                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospy(i), &
3985                                        'units', 'unknown' )
3986                CALL handle_netcdf_error( 'netcdf', 290 )
3987
3988                var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
3989
3990             ENDIF
3991
3992             i = i + 1
3993
3994          ENDDO
3995
3996!
3997!--       Write the list of variables as global attribute (this is used by
3998!--       restart runs)
3999          nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', var_list )
4000          CALL handle_netcdf_error( 'netcdf', 291 )
4001
4002!
4003!--       Leave NetCDF define mode
4004          nc_stat = NF90_ENDDEF( id_set_sp )
4005          CALL handle_netcdf_error( 'netcdf', 292 )
4006
4007!
4008!--       Write axis data: zu_sp, zw_sp, k_x, k_y
4009          ALLOCATE( netcdf_data(1:ns) )
4010
4011!
4012!--       Write zu data
4013          netcdf_data(1:ns) = zu( comp_spectra_level(1:ns) )
4014          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zu_sp, netcdf_data, &
4015                                  start = (/ 1 /), count = (/ ns /) )
4016          CALL handle_netcdf_error( 'netcdf', 293 )
4017
4018!
4019!--       Write zw data
4020          netcdf_data(1:ns) = zw( comp_spectra_level(1:ns) )
4021          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zw_sp, netcdf_data, &
4022                                  start = (/ 1 /), count = (/ ns /) )
4023          CALL handle_netcdf_error( 'netcdf', 294 )
4024
4025          DEALLOCATE( netcdf_data )
4026
4027!
4028!--       Write data for x and y axis (wavenumbers)
4029          ALLOCATE( netcdf_data(nx/2) )
4030          DO  i = 1, nx/2
4031             netcdf_data(i) = 2.0 * pi * i / ( dx * ( nx + 1 ) )
4032          ENDDO
4033
4034          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_x_sp, netcdf_data, &
4035                                  start = (/ 1 /), count = (/ nx/2 /) )
4036          CALL handle_netcdf_error( 'netcdf', 295 )
4037
4038          DEALLOCATE( netcdf_data )
4039
4040          ALLOCATE( netcdf_data(ny/2) )
4041          DO  i = 1, ny/2
4042             netcdf_data(i) = 2.0 * pi * i / ( dy * ( ny + 1 ) )
4043          ENDDO
4044
4045          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_y_sp, netcdf_data, &
4046                                  start = (/ 1 /), count = (/ ny/2 /) )
4047          CALL handle_netcdf_error( 'netcdf', 296 )
4048
4049          DEALLOCATE( netcdf_data )
4050
4051
4052       CASE ( 'sp_ext' )
4053
4054!
4055!--       Get the list of variables and compare with the actual run.
4056!--       First var_list_old has to be reset, since GET_ATT does not assign
4057!--       trailing blanks.
4058          var_list_old = ' '
4059          nc_stat = NF90_GET_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', &
4060                                  var_list_old )
4061          CALL handle_netcdf_error( 'netcdf', 297 )
4062          var_list = ';'
4063          i = 1
4064          DO WHILE ( data_output_sp(i) /= ' '  .AND.  i <= 10 )
4065
4066             IF ( INDEX( spectra_direction(i), 'x' ) /= 0 )  THEN
4067                netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
4068                var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
4069             ENDIF
4070
4071             IF ( INDEX( spectra_direction(i), 'y' ) /= 0 )  THEN
4072                netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
4073                var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
4074             ENDIF
4075
4076             i = i + 1
4077
4078          ENDDO
4079
4080          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
4081             message_string = 'NetCDF file for spectra  ' //                 &
4082                              'from previous run found,' //                  &
4083                              '& but this file cannot be extended due to' // &
4084                              ' variable mismatch.' //                       &
4085                              '&New file is created instead.'
4086             CALL message( 'define_netcdf_header', 'PA0260', 0, 1, 0, 6, 0 )
4087             extend = .FALSE.
4088             RETURN
4089          ENDIF
4090
4091!
4092!--       Determine the number of current vertical levels for which spectra
4093!--       shall be output
4094          ns = 1
4095          DO WHILE ( comp_spectra_level(ns) /= 999999  .AND.  ns <= 100 )
4096             ns = ns + 1
4097          ENDDO
4098          ns = ns - 1
4099
4100!
4101!--       Get and compare the number of vertical levels
4102          nc_stat = NF90_INQ_VARID( id_set_sp, 'zu_sp', id_var_zu_sp )
4103          CALL handle_netcdf_error( 'netcdf', 298 )
4104
4105          nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_zu_sp, &
4106                                           dimids = id_dim_zu_sp_old )
4107          CALL handle_netcdf_error( 'netcdf', 299 )
4108          id_dim_zu_sp = id_dim_zu_sp_old(1)
4109
4110          nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_zu_sp, &
4111                                            len = ns_old )
4112          CALL handle_netcdf_error( 'netcdf', 300 )
4113
4114          IF ( ns /= ns_old )  THEN
4115             message_string = 'NetCDF file for spectra ' //                 &
4116                              ' from previous run found,' //                &
4117                              '&but this file cannot be extended due to' // &
4118                              ' mismatch in number of' //                   &
4119                              '&vertical levels.' //                        &
4120                              '&New file is created instead.'
4121             CALL message( 'define_netcdf_header', 'PA0261', 0, 1, 0, 6, 0 )
4122             extend = .FALSE.
4123             RETURN
4124          ENDIF
4125
4126!
4127!--       Get and compare the heights of the cross sections
4128          ALLOCATE( netcdf_data(1:ns_old) )
4129
4130          nc_stat = NF90_GET_VAR( id_set_sp, id_var_zu_sp, netcdf_data )
4131          CALL handle_netcdf_error( 'netcdf', 301 )
4132
4133          DO  i = 1, ns
4134             IF ( zu(comp_spectra_level(i)) /= netcdf_data(i) )  THEN
4135                message_string = 'NetCDF file for spectra ' //                 &
4136                                 ' from previous run found,' //                &
4137                                 '&but this file cannot be extended due to' // &
4138                                 ' mismatch in heights of' //                  &
4139                                 '&vertical levels.' //                        &
4140                                 '&New file is created instead.'
4141                CALL message( 'define_netcdf_header', 'PA0262', 0, 1, 0, 6, 0 )
4142                extend = .FALSE.
4143                RETURN
4144             ENDIF
4145          ENDDO
4146
4147          DEALLOCATE( netcdf_data )
4148
4149!
4150!--       Get the id of the time coordinate (unlimited coordinate) and its
4151!--       last index on the file. The next time level is plsp..count+1.
4152!--       The current time must be larger than the last output time
4153!--       on the file.
4154          nc_stat = NF90_INQ_VARID( id_set_sp, 'time', id_var_time_sp )
4155          CALL handle_netcdf_error( 'netcdf', 302 )
4156
4157          nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_time_sp, &
4158                                           dimids = id_dim_time_old )
4159          CALL handle_netcdf_error( 'netcdf', 303 )
4160          id_dim_time_sp = id_dim_time_old(1)
4161
4162          nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_time_sp, &
4163                                            len = dosp_time_count )
4164          CALL handle_netcdf_error( 'netcdf', 304 )
4165
4166          nc_stat = NF90_GET_VAR( id_set_sp, id_var_time_sp,        &
4167                                  last_time_coordinate,             &
4168                                  start = (/ dosp_time_count /), &
4169                                  count = (/ 1 /) )
4170          CALL handle_netcdf_error( 'netcdf', 305 )
4171
4172          IF ( last_time_coordinate(1) >= simulated_time )  THEN
4173             message_string = 'NetCDF file for spectra ' //                 &
4174                              'from previous run found,' //                 &
4175                              '&but this file cannot be extended becaus' // &
4176                              'e the current output time' //                & 
4177                              '&is less or equal than the last output t' // &
4178                              'ime on this file.' //                        &
4179                              '&New file is created instead.'
4180             CALL message( 'define_netcdf_header', 'PA0263', 0, 1, 0, 6, 0 )
4181             dosp_time_count = 0
4182             extend = .FALSE.
4183             RETURN
4184          ENDIF
4185
4186!
4187!--       Dataset seems to be extendable.
4188!--       Now get the variable ids.
4189          i = 1
4190          DO WHILE ( data_output_sp(i) /= ' '  .AND.  i <= 10 )
4191
4192             IF ( INDEX( spectra_direction(i), 'x' ) /= 0 )  THEN
4193                netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
4194                nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
4195                                          id_var_dospx(i) )
4196                CALL handle_netcdf_error( 'netcdf', 306 )
4197             ENDIF
4198
4199             IF ( INDEX( spectra_direction(i), 'y' ) /= 0 )  THEN
4200                netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
4201                nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
4202                                          id_var_dospy(i) )
4203                CALL handle_netcdf_error( 'netcdf', 307 )
4204             ENDIF
4205
4206             i = i + 1
4207
4208          ENDDO
4209
4210!
4211!--       Update the title attribute on file
4212!--       In order to avoid 'data mode' errors if updated attributes are larger
4213!--       than their original size, NF90_PUT_ATT is called in 'define mode'
4214!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
4215!--       performance loss due to data copying; an alternative strategy would be
4216!--       to ensure equal attribute size in a job chain. Maybe revise later.
4217          nc_stat = NF90_REDEF( id_set_sp )
4218          CALL handle_netcdf_error( 'netcdf', 441 )
4219          IF ( averaging_interval_sp /= 0.0 )  THEN
4220             WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
4221                                                           averaging_interval_sp
4222             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title',  &
4223                                     TRIM( run_description_header ) // &
4224                                     TRIM( time_average_text ) )
4225             CALL handle_netcdf_error( 'netcdf', 308 )
4226
4227             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval_sp
4228             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'time_avg', &
4229                                     TRIM( time_average_text ) )
4230          ELSE
4231             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
4232                                     TRIM( run_description_header ) )
4233          ENDIF
4234          CALL handle_netcdf_error( 'netcdf', 309 )
4235          nc_stat = NF90_ENDDEF( id_set_sp )
4236          CALL handle_netcdf_error( 'netcdf', 442 )
4237          message_string = 'NetCDF file for spectra ' //     & 
4238                           'from previous run found.' //     &
4239                           '&This file will be extended.'
4240          CALL message( 'define_netcdf_header', 'PA0264', 0, 0, 0, 6, 0 )
4241
4242
4243       CASE ( 'pt_new' )
4244
4245!
4246!--       Define some global attributes of the dataset
4247          nc_stat = NF90_PUT_ATT( id_set_prt, NF90_GLOBAL, 'title', &
4248                                  TRIM( run_description_header ) )
4249          CALL handle_netcdf_error( 'netcdf', 310 )
4250
4251!
4252!--       Define time coordinate for particles (unlimited dimension)
4253          nc_stat = NF90_DEF_DIM( id_set_prt, 'time', NF90_UNLIMITED, &
4254                                  id_dim_time_prt )
4255          CALL handle_netcdf_error( 'netcdf', 311 )
4256
4257          nc_stat = NF90_DEF_VAR( id_set_prt, 'time', NF90_DOUBLE, &
4258                                  id_dim_time_prt, id_var_time_prt )
4259          CALL handle_netcdf_error( 'netcdf', 312 )
4260
4261          nc_stat = NF90_PUT_ATT( id_set_prt, id_var_time_prt, 'units', &
4262                                  'seconds' )
4263          CALL handle_netcdf_error( 'netcdf', 313 )
4264
4265!
4266!--       Define particle coordinate (maximum particle number)
4267          IF ( netcdf_data_format < 3 )  THEN
4268             nc_stat = NF90_DEF_DIM( id_set_prt, 'particle_number', &
4269                                     maximum_number_of_particles, id_dim_prtnum)
4270          ELSE
4271!
4272!--          NetCDF4 allows more than one unlimited dimension
4273             nc_stat = NF90_DEF_DIM( id_set_prt, 'particle_number', &
4274                                     NF90_UNLIMITED, id_dim_prtnum)
4275          ENDIF
4276          CALL handle_netcdf_error( 'netcdf', 314 )
4277
4278          nc_stat = NF90_DEF_VAR( id_set_prt, 'particle_number', NF90_DOUBLE, &
4279                                  id_dim_prtnum, id_var_prtnum )
4280          CALL handle_netcdf_error( 'netcdf', 315 )
4281
4282          nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prtnum, 'units', &
4283                                  'particle number' )
4284          CALL handle_netcdf_error( 'netcdf', 316 )
4285
4286!
4287!--       Define variable which contains the real number of particles in use
4288          nc_stat = NF90_DEF_VAR( id_set_prt, 'real_num_of_prt', NF90_INT, &
4289                                  id_dim_time_prt, id_var_rnop_prt )
4290          CALL handle_netcdf_error( 'netcdf', 317 )
4291
4292          nc_stat = NF90_PUT_ATT( id_set_prt, id_var_rnop_prt, 'units', &
4293                                  'particle number' )
4294          CALL handle_netcdf_error( 'netcdf', 318 )
4295
4296!
4297!--       Define the variables
4298          DO  i = 1, 17
4299
4300             nc_stat = NF90_DEF_VAR( id_set_prt, prt_var_names(i),         &
4301                                     nc_precision(8),                      &
4302                                     (/ id_dim_prtnum, id_dim_time_prt /), &
4303                                     id_var_prt(i) )
4304             CALL handle_netcdf_error( 'netcdf', 319 )
4305
4306             nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prt(i), &
4307                                     'long_name', TRIM( prt_var_names(i) ) )
4308             CALL handle_netcdf_error( 'netcdf', 320 )
4309
4310             nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prt(i), &
4311                                     'units', TRIM( prt_var_units(i) ) )
4312             CALL handle_netcdf_error( 'netcdf', 321 )
4313
4314          ENDDO
4315
4316!
4317!--       Leave NetCDF define mode
4318          nc_stat = NF90_ENDDEF( id_set_prt )
4319          CALL handle_netcdf_error( 'netcdf', 322 )
4320
4321
4322       CASE ( 'pt_ext' )
4323
4324!
4325!--       Get the id of the time coordinate (unlimited coordinate) and its
4326!--       last index on the file. The next time level is prt..count+1.
4327!--       The current time must be larger than the last output time
4328!--       on the file.
4329          nc_stat = NF90_INQ_VARID( id_set_prt, 'time', id_var_time_prt )
4330          CALL handle_netcdf_error( 'netcdf', 323 )
4331
4332          nc_stat = NF90_INQUIRE_VARIABLE( id_set_prt, id_var_time_prt, &
4333                                           dimids = id_dim_time_old )
4334          CALL handle_netcdf_error( 'netcdf', 324 )
4335          id_dim_time_prt = id_dim_time_old(1)
4336
4337          nc_stat = NF90_INQUIRE_DIMENSION( id_set_prt, id_dim_time_prt, &
4338                                            len = prt_time_count )
4339          CALL handle_netcdf_error( 'netcdf', 325 )
4340
4341          nc_stat = NF90_GET_VAR( id_set_prt, id_var_time_prt,  &
4342                                  last_time_coordinate,         &
4343                                  start = (/ prt_time_count /), &
4344                                  count = (/ 1 /) )
4345          CALL handle_netcdf_error( 'netcdf', 326 )
4346
4347          IF ( last_time_coordinate(1) >= simulated_time )  THEN
4348             message_string = 'NetCDF file for particles ' //               &
4349                              'from previous run found,' //                 &
4350                              '&but this file cannot be extended becaus' // &
4351                              'e the current output time' //                &
4352                              '&is less or equal than the last output t' // &
4353                              'ime on this file.' //                        &
4354                              '&New file is created instead.'
4355             CALL message( 'define_netcdf_header', 'PA0265', 0, 1, 0, 6, 0 )
4356             prt_time_count = 0
4357             extend = .FALSE.
4358             RETURN
4359          ENDIF
4360
4361!
4362!--       Dataset seems to be extendable.
4363!--       Now get the variable ids.
4364          nc_stat = NF90_INQ_VARID( id_set_prt, 'real_num_of_prt', &
4365                                    id_var_rnop_prt )
4366          CALL handle_netcdf_error( 'netcdf', 327 )
4367
4368          DO  i = 1, 17
4369
4370             nc_stat = NF90_INQ_VARID( id_set_prt, prt_var_names(i), &
4371                                       id_var_prt(i) )
4372             CALL handle_netcdf_error( 'netcdf', 328 )
4373
4374          ENDDO
4375
4376          message_string = 'NetCDF file for particles ' // & 
4377                           'from previous run found.' //   &
4378                           '&This file will be extended.'
4379          CALL message( 'define_netcdf_header', 'PA0266', 0, 0, 0, 6, 0 )
4380       
4381
4382
4383       CASE ( 'ps_new' )
4384
4385!
4386!--       Define some global attributes of the dataset
4387          nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
4388                                  TRIM( run_description_header ) )
4389          CALL handle_netcdf_error( 'netcdf', 396 )
4390
4391!
4392!--       Define time coordinate for particle time series (unlimited dimension)
4393          nc_stat = NF90_DEF_DIM( id_set_pts, 'time', NF90_UNLIMITED, &
4394                                  id_dim_time_pts )
4395          CALL handle_netcdf_error( 'netcdf', 397 )
4396
4397          nc_stat = NF90_DEF_VAR( id_set_pts, 'time', NF90_DOUBLE, &
4398                                  id_dim_time_pts, id_var_time_pts )
4399          CALL handle_netcdf_error( 'netcdf', 398 )
4400
4401          nc_stat = NF90_PUT_ATT( id_set_pts, id_var_time_pts, 'units', &
4402                                  'seconds')
4403          CALL handle_netcdf_error( 'netcdf', 399 )
4404
4405!
4406!--       Define the variables. If more than one particle group is defined,
4407!--       define seperate variables for each group
4408          var_list = ';'
4409          DO  i = 1, dopts_num
4410
4411             DO  j = 0, number_of_particle_groups
4412
4413                IF ( j == 0 )  THEN
4414                   suffix1 = ''
4415                ELSE
4416                   WRITE ( suffix1, '(''_'',I2.2)' )  j
4417                ENDIF
4418
4419                nc_stat = NF90_DEF_VAR( id_set_pts,                         &
4420                                        TRIM( dopts_label(i) ) // suffix1,  &
4421                                        nc_precision(6),                    &
4422                                        (/ id_dim_time_pts /),              &
4423                                        id_var_dopts(i,j) )
4424                CALL handle_netcdf_error( 'netcdf', 400 )
4425
4426                IF ( j == 0 )  THEN
4427                   nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
4428                                           'long_name', &
4429                                           TRIM( dopts_label(i) ) )
4430                ELSE
4431                   nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
4432                                           'long_name', &
4433                                           TRIM( dopts_label(i) ) // ' PG ' // &
4434                                           suffix1(2:3) )
4435                ENDIF
4436                CALL handle_netcdf_error( 'netcdf', 401 )
4437
4438                nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
4439                                        'units', TRIM( dopts_unit(i) ) )
4440                CALL handle_netcdf_error( 'netcdf', 402 )
4441
4442                var_list = TRIM( var_list ) // TRIM( dopts_label(i) ) // &
4443                           suffix1 // ';'
4444
4445                IF ( number_of_particle_groups == 1 )  EXIT
4446
4447             ENDDO
4448
4449          ENDDO
4450
4451!
4452!--       Write the list of variables as global attribute (this is used by
4453!--       restart runs)
4454          nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
4455                                  var_list )
4456          CALL handle_netcdf_error( 'netcdf', 403 )
4457
4458
4459!
4460!--       Leave NetCDF define mode
4461          nc_stat = NF90_ENDDEF( id_set_pts )
4462          CALL handle_netcdf_error( 'netcdf', 404 )
4463
4464
4465       CASE ( 'ps_ext' )
4466
4467!
4468!--       Get the list of variables and compare with the actual run.
4469!--       First var_list_old has to be reset, since GET_ATT does not assign
4470!--       trailing blanks.
4471          var_list_old = ' '
4472          nc_stat = NF90_GET_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
4473                                  var_list_old )
4474          CALL handle_netcdf_error( 'netcdf', 405 )
4475
4476          var_list = ';'
4477          i = 1
4478          DO  i = 1, dopts_num
4479
4480             DO  j = 0, number_of_particle_groups
4481
4482                IF ( j == 0 )  THEN
4483                   suffix1 = ''
4484                ELSE
4485                   WRITE ( suffix1, '(''_'',I2.2)' )  j
4486                ENDIF
4487
4488                var_list = TRIM( var_list ) // TRIM( dopts_label(i) ) // &
4489                           suffix1 // ';'
4490
4491                IF ( number_of_particle_groups == 1 )  EXIT
4492
4493             ENDDO
4494
4495          ENDDO
4496
4497          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
4498             message_string = 'NetCDF file for particle time series ' //     &
4499                              'from previous run found,' //                  &
4500                              '& but this file cannot be extended due to' // &
4501                              ' variable mismatch.' //                       &
4502                              '&New file is created instead.'
4503             CALL message( 'define_netcdf_header', 'PA0267', 0, 1, 0, 6, 0 )
4504             extend = .FALSE.
4505             RETURN
4506          ENDIF
4507
4508!
4509!--       Get the id of the time coordinate (unlimited coordinate) and its
4510!--       last index on the file. The next time level is dots..count+1.
4511!--       The current time must be larger than the last output time
4512!--       on the file.
4513          nc_stat = NF90_INQ_VARID( id_set_pts, 'time', id_var_time_pts )
4514          CALL handle_netcdf_error( 'netcdf', 406 )
4515
4516          nc_stat = NF90_INQUIRE_VARIABLE( id_set_pts, id_var_time_pts, &
4517                                           dimids = id_dim_time_old )
4518          CALL handle_netcdf_error( 'netcdf', 407 )
4519          id_dim_time_pts = id_dim_time_old(1)
4520
4521          nc_stat = NF90_INQUIRE_DIMENSION( id_set_pts, id_dim_time_pts, &
4522                                            len = dopts_time_count )
4523          CALL handle_netcdf_error( 'netcdf', 408 )
4524
4525          nc_stat = NF90_GET_VAR( id_set_pts, id_var_time_pts,    &
4526                                  last_time_coordinate,           &
4527                                  start = (/ dopts_time_count /), &
4528                                  count = (/ 1 /) )
4529          CALL handle_netcdf_error( 'netcdf', 409 )
4530
4531          IF ( last_time_coordinate(1) >= simulated_time )  THEN
4532             message_string = 'NetCDF file for particle time series ' //    &
4533                              'from previous run found,' //                 &
4534                              '&but this file cannot be extended becaus' // &
4535                              'e the current output time' //                &
4536                              '&is less or equal than the last output t' // &
4537                              'ime on this file.' //                        &
4538                              '&New file is created instead.'
4539             CALL message( 'define_netcdf_header', 'PA0268', 0, 1, 0, 6, 0 )
4540             dopts_time_count = 0
4541             extend = .FALSE.
4542             RETURN
4543          ENDIF
4544
4545!
4546!--       Dataset seems to be extendable.
4547!--       Now get the variable ids
4548          i = 1
4549          DO  i = 1, dopts_num
4550 
4551             DO  j = 0, number_of_particle_groups
4552
4553                IF ( j == 0 )  THEN
4554                   suffix1 = ''
4555                ELSE
4556                   WRITE ( suffix1, '(''_'',I2.2)' )  j
4557                ENDIF
4558
4559                netcdf_var_name = TRIM( dopts_label(i) ) // suffix1
4560
4561                nc_stat = NF90_INQ_VARID( id_set_pts, netcdf_var_name, &
4562                                          id_var_dopts(i,j) )
4563                CALL handle_netcdf_error( 'netcdf', 410 )
4564
4565                IF ( number_of_particle_groups == 1 )  EXIT
4566
4567             ENDDO
4568
4569          ENDDO
4570
4571!
4572!--       Update the title attribute on file
4573!--       In order to avoid 'data mode' errors if updated attributes are larger
4574!--       than their original size, NF90_PUT_ATT is called in 'define mode'
4575!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
4576!--       performance loss due to data copying; an alternative strategy would be
4577!--       to ensure equal attribute size in a job chain. Maybe revise later.
4578          nc_stat = NF90_REDEF( id_set_pts )
4579          CALL handle_netcdf_error( 'netcdf', 443 )
4580          nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
4581                                  TRIM( run_description_header ) )
4582          CALL handle_netcdf_error( 'netcdf', 411 )
4583          nc_stat = NF90_ENDDEF( id_set_pts )
4584          CALL handle_netcdf_error( 'netcdf', 444 )
4585          message_string = 'NetCDF file for particle time series ' // & 
4586                           'from previous run found.' //              &
4587                           '&This file will be extended.'
4588          CALL message( 'define_netcdf_header', 'PA0269', 0, 0, 0, 6, 0 )
4589
4590
4591       CASE DEFAULT
4592
4593          message_string = 'mode "' // TRIM( mode) // '" not supported'
4594          CALL message( 'define_netcdf_header', 'PA0270', 0, 0, 0, 6, 0 )
4595
4596    END SELECT
4597
4598#endif
4599 END SUBROUTINE define_netcdf_header
4600
4601
4602SUBROUTINE handle_netcdf_error( routine_name, errno )
4603#if defined( __netcdf )
4604
4605!------------------------------------------------------------------------------!
4606!
4607! Description:
4608! ------------
4609! Prints out a text message corresponding to the current status.
4610!------------------------------------------------------------------------------!
4611
4612    USE control_parameters
4613    USE netcdf
4614    USE netcdf_control
4615    USE pegrid
4616
4617    IMPLICIT NONE
4618
4619    CHARACTER(LEN=6) ::  message_identifier
4620    CHARACTER(LEN=*) ::  routine_name
4621
4622    INTEGER ::  errno
4623
4624    IF ( nc_stat /= NF90_NOERR )  THEN
4625
4626       WRITE( message_identifier, '(''NC'',I4.4)' )  errno
4627       message_string = TRIM( NF90_STRERROR( nc_stat ) )
4628
4629       CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
4630
4631    ENDIF
4632
4633#endif
4634 END SUBROUTINE handle_netcdf_error
Note: See TracBrowser for help on using the repository browser.