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

Last change on this file since 1036 was 1036, checked in by raasch, 11 years ago

code has been put under the GNU General Public License (v3)

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