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

Last change on this file since 1054 was 1054, checked in by hoffmann, 11 years ago

last commit documented

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