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

Last change on this file since 1585 was 1585, checked in by maronga, 6 years ago

Added support for RRTMG radiation code

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