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

Last change on this file since 1093 was 1093, checked in by raasch, 9 years ago

last commit documented

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