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

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

last commit documented

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