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

Last change on this file since 1552 was 1552, checked in by maronga, 7 years ago

last commit documented

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