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

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

land surface model released

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