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

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

New:
---
Output in NetCDF4-format. New d3par-parameter netcdf_data_format.

(check_open, check_parameters, close_file, data_output_2d, data_output_3d, header, modules, netcdf, parin)

Modules to be loaded for compilation (mbuild) or job execution (mrun)
can be given in the configuration file using variable modules. Example:

%modules ifort/11.0.069:netcdf lcsgih parallel

This method replaces the (undocumented) mpilib-variable.

WARNING: All fixed settings of modules in the scripts mbuild, mrun, and subjob
have been removed! Please set the modules variable appropriately in your
configuration file. (mbuild, mrun, subjob)

Changed:


Parameters netcdf_64bit and netcdf_64bit_3d have been removed. Use
netcdf_data_format = 2 for choosing the classic 64bit-offset format (this is
the default). The offset-format can not be set independently for the
3d-output-data any more.

Parameters netcdf_format_mask, netcdf_format_mask_av, and variables
nc_format_mask, format_parallel_io removed. They are replaced by the new
parameter netcdf_data_format. (check_open, close_file,
data_output_mask, header, init_masks, modules, parin)

Errors:


bugfix in trunk/UTIL/Makefile: forgot to compile for interpret_config

Bugfix: timeseries data have to be collected by PE0 (user_statistics)

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