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

Last change on this file since 953 was 952, checked in by hoffmann, 12 years ago

last commit documented

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