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

Last change on this file since 1207 was 1207, checked in by witha, 11 years ago

last commit documented

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