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

Last change on this file since 1309 was 1309, checked in by fricke, 7 years ago

last commit documented

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