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

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

various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings

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