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

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

NetCDF4 support for particle data; special character allowed for NetCDF variable names

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