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

Last change on this file since 1683 was 1683, checked in by knoop, 6 years ago

last commit documented

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