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

Last change on this file since 1321 was 1321, checked in by raasch, 7 years ago

last commit documented

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