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

Last change on this file since 964 was 964, checked in by raasch, 12 years ago

old profil-parameters (cross_xtext, cross_normalized_x, etc. ) and respective code removed
(check_open, check_parameters, close_file, data_output_profiles, data_output_spectra, header, modules, parin)

reformatting (netcdf)

append feature removed from unit 14 (check_open)

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