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

Last change on this file since 960 was 960, checked in by hoffmann, 9 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 183.9 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!
10!
11! Former revisions:
12! -----------------
13! $Id: netcdf.f90 960 2012-07-24 13:21:34Z hoffmann $
14!
15! 959 2012-07-24 13:13:41Z hoffmann
16! Bugfix in cross_profiles. It is not possible 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, 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 columns and rows of cross_profiles to netcdf header.
3150!--       This information can be used by palmplot.
3151
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          k = 1
3167
3168!
3169!--       Each profile defined in cross_profiles is written to a array
3170!--       (cross_profiles_char). The number of its cross is assigned in a
3171!--       second array (cross_profiles_numb).
3172
3173          DO  i = 1, crmax
3174            IF ( TRIM( cross_profiles_adj(i) ) == ' ' )  EXIT
3175            delim_old = 0
3176            DO   j = 1, crmax
3177              delim = INDEX( cross_profiles_adj(i)(delim_old+1:), ' ' )
3178              IF (delim .EQ. 1)  EXIT
3179              cross_profiles_char( MIN( crmax, k ) ) = cross_profiles_adj(i)  &
3180                                       (delim_old+1:delim_old+delim-1)
3181              cross_profiles_numb( MIN( crmax, k ) ) = i
3182              k = k+1
3183              cross_profiles_maxi  = i
3184              delim_old = delim_old + delim
3185            ENDDO
3186          ENDDO
3187
3188          cross_profiles_count = MIN( crmax, k-1 )
3189         
3190!
3191!--       Check if all profiles defined in cross_profiles are defined in
3192!--       data_output_pr. If not, they will be skipped.
3193
3194          cross_profiles_log = .FALSE.
3195          message_string = ' '
3196
3197          DO  i = 1, cross_profiles_count
3198            DO  j = 1, dopr_n
3199              IF ( TRIM(cross_profiles_char(i) ) == TRIM(data_output_pr(j)) ) &
3200                 EXIT
3201              IF ( j == dopr_n)  THEN
3202                cross_profiles_numb(i) = 999999
3203                cross_profiles_log = .TRUE.
3204                message_string = TRIM(message_string) // ', ' //              &
3205                                 TRIM(cross_profiles_char(i))
3206              ENDIF
3207            ENDDO
3208          ENDDO
3209
3210          IF (cross_profiles_log)  THEN
3211            message_string = TRIM(message_string(2:)) // ' is/are not' //     &
3212                             ' defined in data_output_pr.'
3213            CALL message( 'define_netcdf_header', 'PA0352', 0, 0, 0, 6, 0 )
3214          ENDIF
3215
3216!
3217!--       Check if all profiles defined in data_output_pr are defined in
3218!--       cross_profiles. If not, they will be added to cross_profiles.
3219
3220          cross_profiles_log = .FALSE.
3221          message_string = ' '
3222
3223          DO  i = 1, dopr_n
3224            DO  j = 1, cross_profiles_count
3225              IF ( TRIM( cross_profiles_char(j) ) == TRIM( data_output_pr(i)) )&
3226                 EXIT
3227              IF ( j == cross_profiles_count )  THEN
3228                cross_profiles_count = cross_profiles_count + 1
3229                cross_profiles_maxi  = cross_profiles_maxi  + 1
3230                cross_profiles_char( MIN( crmax, cross_profiles_count ) ) =    &
3231                                                      TRIM( data_output_pr(i) )
3232                cross_profiles_numb( MIN( crmax, cross_profiles_count ) ) =    &
3233                                                      cross_profiles_maxi
3234                cross_profiles_log = .TRUE.
3235                message_string = TRIM( message_string ) // ', '&
3236                                         // TRIM( data_output_pr(i) )
3237              ENDIF
3238            ENDDO
3239          ENDDO
3240
3241          IF ( cross_profiles_log )  THEN
3242            message_string = TRIM(message_string(2:)) //               &
3243                             ' has/have been added to cross_profiles.'
3244            CALL message( 'define_netcdf_header', 'PA0353', 0, 0, 0, 6, 0 )
3245          ENDIF
3246
3247          IF ( cross_profiles_count .ge. crmax )  THEN
3248            message_string = 'It is not possible to arrange more than '&
3249                             // '100 profiles with cross_profiles.'
3250            CALL message( 'define_netcdf_header', 'PA0354', 0, 0, 0, 6, 0 )
3251          ENDIF
3252
3253!
3254!--       Writing cross_profiles to netcdf header. This information can be
3255!--       used by palmplot. Each profile is separated by ",", each cross is
3256!--       separated by ";".
3257
3258          char_cross_profiles = ';'
3259          id_last = 1
3260
3261          DO  i = 1, cross_profiles_count
3262            IF ( cross_profiles_numb(i) /= 999999 )  THEN
3263              IF ( TRIM( char_cross_profiles ) == ';' )  THEN
3264               char_cross_profiles = TRIM( char_cross_profiles ) // &
3265                                     TRIM( cross_profiles_char(i) )
3266              ELSEIF ( id_last == cross_profiles_numb(i) )  THEN
3267                char_cross_profiles = TRIM( char_cross_profiles ) // &
3268                                      ',' // TRIM( cross_profiles_char(i) )
3269              ELSE
3270                char_cross_profiles = TRIM( char_cross_profiles ) // &
3271                                      ';' // TRIM( cross_profiles_char(i) )
3272              ENDIF
3273              id_last = cross_profiles_numb(i)
3274            ENDIF
3275          ENDDO
3276         
3277          char_cross_profiles = TRIM( char_cross_profiles ) // ';'
3278
3279          nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'cross_profiles',   &
3280                                  TRIM( char_cross_profiles ) )
3281          CALL handle_netcdf_error( 'netcdf', 521 )
3282
3283!
3284!--       Define time coordinate for profiles (unlimited dimension)
3285          nc_stat = NF90_DEF_DIM( id_set_pr, 'time', NF90_UNLIMITED, &
3286                                  id_dim_time_pr )
3287          CALL handle_netcdf_error( 'netcdf', 220 )
3288
3289          nc_stat = NF90_DEF_VAR( id_set_pr, 'time', NF90_DOUBLE, &
3290                                  id_dim_time_pr, id_var_time_pr )
3291          CALL handle_netcdf_error( 'netcdf', 221 )
3292
3293          nc_stat = NF90_PUT_ATT( id_set_pr, id_var_time_pr, 'units', 'seconds')
3294          CALL handle_netcdf_error( 'netcdf', 222 )
3295
3296!
3297!--       Define the variables
3298          var_list = ';'
3299          DO  i = 1, dopr_n
3300
3301             IF ( statistic_regions == 0 )  THEN
3302
3303!
3304!--             Define the z-axes (each variable gets its own z-axis)
3305                nc_stat = NF90_DEF_DIM( id_set_pr,                        &
3306                                        'z' // TRIM( data_output_pr(i) ), &
3307                                        nzt+2-nzb, id_dim_z_pr(i,0) )
3308                CALL handle_netcdf_error( 'netcdf', 223 )
3309
3310                nc_stat = NF90_DEF_VAR( id_set_pr,                        &
3311                                        'z' // TRIM( data_output_pr(i) ), &
3312                                        NF90_DOUBLE, id_dim_z_pr(i,0),    &
3313                                        id_var_z_pr(i,0) )
3314                CALL handle_netcdf_error( 'netcdf', 224 )
3315
3316                nc_stat = NF90_PUT_ATT( id_set_pr, id_var_z_pr(i,0), 'units', &
3317                                        'meters' )
3318                CALL handle_netcdf_error( 'netcdf', 225 )
3319!
3320!--             Define the variable
3321                nc_stat = NF90_DEF_VAR( id_set_pr, data_output_pr(i),         &
3322                                        nc_precision(5), (/ id_dim_z_pr(i,0), &
3323                                        id_dim_time_pr /), id_var_dopr(i,0) )
3324                CALL handle_netcdf_error( 'netcdf', 226 )
3325
3326                nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,0), &
3327                                        'long_name', TRIM( data_output_pr(i) ) )
3328                CALL handle_netcdf_error( 'netcdf', 227 )
3329
3330                nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,0), &
3331                                        'units', TRIM( dopr_unit(i) ) )
3332                CALL handle_netcdf_error( 'netcdf', 228 )
3333
3334                var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) //  ';'
3335
3336             ELSE
3337!
3338!--             If statistic regions are defined, add suffix _SR+#SR to the
3339!--             names
3340                DO  j = 0, statistic_regions
3341                   WRITE ( suffix, '(''_'',I1)' )  j
3342
3343!
3344!--                Define the z-axes (each variable gets it own z-axis)
3345                   nc_stat = NF90_DEF_DIM( id_set_pr,                          &
3346                                         'z'//TRIM(data_output_pr(i))//suffix, &
3347                                           nzt+2-nzb, id_dim_z_pr(i,j) )
3348                   CALL handle_netcdf_error( 'netcdf', 229 )
3349
3350                   nc_stat = NF90_DEF_VAR( id_set_pr,                          &
3351                                         'z'//TRIM(data_output_pr(i))//suffix, &
3352                                           nc_precision(5), id_dim_z_pr(i,j),  &
3353                                           id_var_z_pr(i,j) )
3354                   CALL handle_netcdf_error( 'netcdf', 230 )
3355
3356                   nc_stat = NF90_PUT_ATT( id_set_pr, id_var_z_pr(i,j), &
3357                                           'units', 'meters' )
3358                   CALL handle_netcdf_error( 'netcdf', 231 )
3359
3360!
3361!--                Define the variable
3362                   nc_stat = NF90_DEF_VAR( id_set_pr,                         &
3363                                           TRIM(data_output_pr(i)) // suffix, &
3364                                           nc_precision(5),                   &
3365                                           (/ id_dim_z_pr(i,j),               &
3366                                           id_dim_time_pr /), id_var_dopr(i,j) )
3367                   CALL handle_netcdf_error( 'netcdf', 232 )
3368
3369                   nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,j),        &
3370                                           'long_name',                        &
3371                                           TRIM( data_output_pr(i) ) // ' SR ' &
3372                                           // suffix(2:2) )
3373                   CALL handle_netcdf_error( 'netcdf', 233 )
3374
3375                   nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,j), &
3376                                           'units', TRIM( dopr_unit(i) ) )
3377                   CALL handle_netcdf_error( 'netcdf', 234 )
3378
3379                   var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) // &
3380                              suffix // ';'
3381
3382                ENDDO
3383
3384             ENDIF
3385
3386          ENDDO
3387
3388!
3389!--       Write the list of variables as global attribute (this is used by
3390!--       restart runs)
3391          nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'VAR_LIST', var_list )
3392          CALL handle_netcdf_error( 'netcdf', 235 )
3393
3394!
3395!--       Define normalization variables (as time series)
3396          DO  i = 1, dopr_norm_num
3397
3398             nc_stat = NF90_DEF_VAR( id_set_pr, 'NORM_' // &
3399                                     TRIM( dopr_norm_names(i) ), &
3400                                     nc_precision(5), (/ id_dim_time_pr /), &
3401                                     id_var_norm_dopr(i) )
3402             CALL handle_netcdf_error( 'netcdf', 236 )
3403
3404             nc_stat = NF90_PUT_ATT( id_set_pr, id_var_norm_dopr(i), &
3405                                     'long_name',                    &
3406                                     TRIM( dopr_norm_longnames(i) ) )
3407             CALL handle_netcdf_error( 'netcdf', 237 )
3408
3409          ENDDO
3410
3411!
3412!--       Leave NetCDF define mode
3413          nc_stat = NF90_ENDDEF( id_set_pr )
3414          CALL handle_netcdf_error( 'netcdf', 238 )
3415
3416!
3417!--       Write z-axes data
3418          DO  i = 1, dopr_n
3419             DO  j = 0, statistic_regions
3420
3421                nc_stat = NF90_PUT_VAR( id_set_pr, id_var_z_pr(i,j),      &
3422                                        hom(nzb:nzt+1,2,dopr_index(i),0), &
3423                                        start = (/ 1 /),                  &
3424                                        count = (/ nzt-nzb+2 /) )
3425                CALL handle_netcdf_error( 'netcdf', 239 )
3426
3427             ENDDO
3428          ENDDO
3429
3430
3431       CASE ( 'pr_ext' )
3432
3433!
3434!--       Get the list of variables and compare with the actual run.
3435!--       First var_list_old has to be reset, since GET_ATT does not assign
3436!--       trailing blanks.
3437          var_list_old = ' '
3438          nc_stat = NF90_GET_ATT( id_set_pr, NF90_GLOBAL, 'VAR_LIST', &
3439                                  var_list_old )
3440          CALL handle_netcdf_error( 'netcdf', 240 )
3441
3442          var_list = ';'
3443          DO  i = 1, dopr_n
3444
3445             IF ( statistic_regions == 0 )  THEN
3446                var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) // ';'
3447             ELSE
3448                DO  j = 0, statistic_regions
3449                   WRITE ( suffix, '(''_'',I1)' )  j
3450                   var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) // &
3451                              suffix // ';'
3452                ENDDO
3453             ENDIF
3454
3455          ENDDO
3456
3457          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
3458             message_string = 'NetCDF file for vertical profiles ' //        &
3459                              'from previous run found,' //                  &
3460                              '& but this file cannot be extended due to' // &
3461                              ' variable mismatch.' //                       &
3462                              '&New file is created instead.'
3463             CALL message( 'define_netcdf_header', 'PA0254', 0, 1, 0, 6, 0 )
3464             extend = .FALSE.
3465             RETURN
3466          ENDIF
3467
3468!
3469!--       Get the id of the time coordinate (unlimited coordinate) and its
3470!--       last index on the file. The next time level is dopr..count+1.
3471!--       The current time must be larger than the last output time
3472!--       on the file.
3473          nc_stat = NF90_INQ_VARID( id_set_pr, 'time', id_var_time_pr )
3474          CALL handle_netcdf_error( 'netcdf', 241 )
3475
3476          nc_stat = NF90_INQUIRE_VARIABLE( id_set_pr, id_var_time_pr, &
3477                                           dimids = id_dim_time_old )
3478          CALL handle_netcdf_error( 'netcdf', 242 )
3479          id_dim_time_pr = id_dim_time_old(1)
3480
3481          nc_stat = NF90_INQUIRE_DIMENSION( id_set_pr, id_dim_time_pr, &
3482                                            len = dopr_time_count )
3483          CALL handle_netcdf_error( 'netcdf', 243 )
3484
3485          nc_stat = NF90_GET_VAR( id_set_pr, id_var_time_pr,        &
3486                                  last_time_coordinate,             &
3487                                  start = (/ dopr_time_count /), &
3488                                  count = (/ 1 /) )
3489          CALL handle_netcdf_error( 'netcdf', 244 )
3490
3491          IF ( last_time_coordinate(1) >= simulated_time )  THEN
3492             message_string = 'NetCDF file for vertical profiles ' //       &
3493                              'from previous run found,' //                 &
3494                              '&but this file cannot be extended becaus' // &
3495                              'e the current output time' //                &
3496                              '&is less or equal than the last output t' // &
3497                              'ime on this file.' //                        &
3498                              '&New file is created instead.'
3499             CALL message( 'define_netcdf_header', 'PA0255', 0, 1, 0, 6, 0 )
3500             dopr_time_count = 0
3501             extend = .FALSE.
3502             RETURN
3503          ENDIF
3504
3505!
3506!--       Dataset seems to be extendable.
3507!--       Now get the variable ids.
3508          i = 1
3509          DO  i = 1, dopr_n
3510 
3511             IF ( statistic_regions == 0 )  THEN
3512                nc_stat = NF90_INQ_VARID( id_set_pr, data_output_pr(i), &
3513                                          id_var_dopr(i,0) )
3514                CALL handle_netcdf_error( 'netcdf', 245 )
3515             ELSE
3516                DO  j = 0, statistic_regions
3517                   WRITE ( suffix, '(''_'',I1)' )  j
3518                   netcdf_var_name = TRIM( data_output_pr(i) ) // suffix
3519                   nc_stat = NF90_INQ_VARID( id_set_pr, netcdf_var_name, &
3520                                             id_var_dopr(i,j) )
3521                   CALL handle_netcdf_error( 'netcdf', 246 )
3522                ENDDO
3523             ENDIF
3524
3525          ENDDO
3526
3527!
3528!--       Get ids of the normalization variables
3529          DO  i = 1, dopr_norm_num
3530             nc_stat = NF90_INQ_VARID( id_set_pr,                             &
3531                                       'NORM_' // TRIM( dopr_norm_names(i) ), &
3532                                       id_var_norm_dopr(i) )
3533             CALL handle_netcdf_error( 'netcdf', 247 )
3534          ENDDO
3535
3536!
3537!--       Update the title attribute on file
3538!--       In order to avoid 'data mode' errors if updated attributes are larger
3539!--       than their original size, NF90_PUT_ATT is called in 'define mode'
3540!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
3541!--       performance loss due to data copying; an alternative strategy would be
3542!--       to ensure equal attribute size in a job chain. Maybe revise later.
3543          IF ( averaging_interval_pr == 0.0 )  THEN
3544             time_average_text = ' '
3545          ELSE
3546             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
3547                                                            averaging_interval_pr
3548          ENDIF
3549          nc_stat = NF90_REDEF( id_set_pr )
3550          CALL handle_netcdf_error( 'netcdf', 437 )
3551          nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'title', &
3552                                  TRIM( run_description_header ) //    &
3553                                  TRIM( time_average_text ) )
3554          CALL handle_netcdf_error( 'netcdf', 248 )
3555
3556          nc_stat = NF90_ENDDEF( id_set_pr )
3557          CALL handle_netcdf_error( 'netcdf', 438 )
3558          message_string = 'NetCDF file for vertical profiles ' // & 
3559                           'from previous run found.' //           &
3560                           '&This file will be extended.'
3561          CALL message( 'define_netcdf_header', 'PA0256', 0, 0, 0, 6, 0 )
3562
3563
3564       CASE ( 'ts_new' )
3565
3566!
3567!--       Define some global attributes of the dataset
3568          nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'title', &
3569                                  TRIM( run_description_header ) )
3570          CALL handle_netcdf_error( 'netcdf', 249 )
3571
3572!
3573!--       Define time coordinate for time series (unlimited dimension)
3574          nc_stat = NF90_DEF_DIM( id_set_ts, 'time', NF90_UNLIMITED, &
3575                                  id_dim_time_ts )
3576          CALL handle_netcdf_error( 'netcdf', 250 )
3577
3578          nc_stat = NF90_DEF_VAR( id_set_ts, 'time', NF90_DOUBLE, &
3579                                  id_dim_time_ts, id_var_time_ts )
3580          CALL handle_netcdf_error( 'netcdf', 251 )
3581
3582          nc_stat = NF90_PUT_ATT( id_set_ts, id_var_time_ts, 'units', 'seconds')
3583          CALL handle_netcdf_error( 'netcdf', 252 )
3584
3585!
3586!--       Define the variables
3587          var_list = ';'
3588          DO  i = 1, dots_num
3589
3590             IF ( statistic_regions == 0 )  THEN
3591
3592                nc_stat = NF90_DEF_VAR( id_set_ts, dots_label(i),              &
3593                                        nc_precision(6), (/ id_dim_time_ts /), &
3594                                        id_var_dots(i,0) )
3595                CALL handle_netcdf_error( 'netcdf', 253 )
3596
3597                nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,0), &
3598                                        'long_name', TRIM( dots_label(i) ) )
3599                CALL handle_netcdf_error( 'netcdf', 254 )
3600
3601                nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,0), &
3602                                        'units', TRIM( dots_unit(i) ) )
3603                CALL handle_netcdf_error( 'netcdf', 255 )
3604
3605                var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // ';'
3606
3607             ELSE
3608!
3609!--             If statistic regions are defined, add suffix _SR+#SR to the
3610!--             names
3611                DO  j = 0, statistic_regions
3612                   WRITE ( suffix, '(''_'',I1)' )  j
3613
3614                   nc_stat = NF90_DEF_VAR( id_set_ts,                         &
3615                                           TRIM( dots_label(i) ) // suffix,   &
3616                                           nc_precision(6),                   &
3617                                           (/ id_dim_time_ts /),              &
3618                                           id_var_dots(i,j) )
3619                   CALL handle_netcdf_error( 'netcdf', 256 )
3620
3621                   nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,j),       &
3622                                           'long_name',                       &
3623                                           TRIM( dots_label(i) ) // ' SR ' // &
3624                                           suffix(2:2) )
3625                   CALL handle_netcdf_error( 'netcdf', 257 )
3626
3627                   nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,j), &
3628                                           'units', TRIM( dots_unit(i) ) )
3629                   CALL handle_netcdf_error( 'netcdf', 347 )
3630
3631                   var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // &
3632                              suffix // ';'
3633
3634                ENDDO
3635
3636             ENDIF
3637
3638          ENDDO
3639
3640!
3641!--       Write the list of variables as global attribute (this is used by
3642!--       restart runs)
3643          nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'VAR_LIST', var_list )
3644          CALL handle_netcdf_error( 'netcdf', 258 )
3645
3646!
3647!--       Leave NetCDF define mode
3648          nc_stat = NF90_ENDDEF( id_set_ts )
3649          CALL handle_netcdf_error( 'netcdf', 259 )
3650
3651
3652       CASE ( 'ts_ext' )
3653
3654!
3655!--       Get the list of variables and compare with the actual run.
3656!--       First var_list_old has to be reset, since GET_ATT does not assign
3657!--       trailing blanks.
3658          var_list_old = ' '
3659          nc_stat = NF90_GET_ATT( id_set_ts, NF90_GLOBAL, 'VAR_LIST', &
3660                                  var_list_old )
3661          CALL handle_netcdf_error( 'netcdf', 260 )
3662
3663          var_list = ';'
3664          i = 1
3665          DO  i = 1, dots_num
3666
3667             IF ( statistic_regions == 0 )  THEN
3668                var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // ';'
3669             ELSE
3670                DO  j = 0, statistic_regions
3671                   WRITE ( suffix, '(''_'',I1)' )  j
3672                   var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // &
3673                              suffix // ';'
3674                ENDDO
3675             ENDIF
3676
3677          ENDDO
3678
3679          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
3680             message_string = 'NetCDF file for time series ' //              &
3681                              'from previous run found,' //                  &
3682                              '& but this file cannot be extended due to' // &
3683                              ' variable mismatch.' //                       &
3684                              '&New file is created instead.'
3685             CALL message( 'define_netcdf_header', 'PA0257', 0, 1, 0, 6, 0 )
3686             extend = .FALSE.
3687             RETURN
3688          ENDIF
3689
3690!
3691!--       Get the id of the time coordinate (unlimited coordinate) and its
3692!--       last index on the file. The next time level is dots..count+1.
3693!--       The current time must be larger than the last output time
3694!--       on the file.
3695          nc_stat = NF90_INQ_VARID( id_set_ts, 'time', id_var_time_ts )
3696          CALL handle_netcdf_error( 'netcdf', 261 )
3697
3698          nc_stat = NF90_INQUIRE_VARIABLE( id_set_ts, id_var_time_ts, &
3699                                           dimids = id_dim_time_old )
3700          CALL handle_netcdf_error( 'netcdf', 262 )
3701          id_dim_time_ts = id_dim_time_old(1)
3702
3703          nc_stat = NF90_INQUIRE_DIMENSION( id_set_ts, id_dim_time_ts, &
3704                                            len = dots_time_count )
3705          CALL handle_netcdf_error( 'netcdf', 263 )
3706
3707          nc_stat = NF90_GET_VAR( id_set_ts, id_var_time_ts,        &
3708                                  last_time_coordinate,             &
3709                                  start = (/ dots_time_count /), &
3710                                  count = (/ 1 /) )
3711          CALL handle_netcdf_error( 'netcdf', 264 )
3712
3713          IF ( last_time_coordinate(1) >= simulated_time )  THEN
3714             message_string = 'NetCDF file for time series ' //             &
3715                              'from previous run found,' //                 &
3716                              '&but this file cannot be extended becaus' // &
3717                              'e the current output time' //                &
3718                              '&is less or equal than the last output t' // &
3719                              'ime on this file.' //                        &
3720                              '&New file is created instead.'
3721             CALL message( 'define_netcdf_header', 'PA0258', 0, 1, 0, 6, 0 )
3722             dots_time_count = 0
3723             extend = .FALSE.
3724             RETURN
3725          ENDIF
3726
3727!
3728!--       Dataset seems to be extendable.
3729!--       Now get the variable ids
3730          i = 1
3731          DO  i = 1, dots_num
3732 
3733             IF ( statistic_regions == 0 )  THEN
3734                nc_stat = NF90_INQ_VARID( id_set_ts, dots_label(i), &
3735                                          id_var_dots(i,0) )
3736                CALL handle_netcdf_error( 'netcdf', 265 )
3737             ELSE
3738                DO  j = 0, statistic_regions
3739                   WRITE ( suffix, '(''_'',I1)' )  j
3740                   netcdf_var_name = TRIM( dots_label(i) ) // suffix
3741                   nc_stat = NF90_INQ_VARID( id_set_ts, netcdf_var_name, &
3742                                             id_var_dots(i,j) )
3743                   CALL handle_netcdf_error( 'netcdf', 266 )
3744                ENDDO
3745             ENDIF
3746
3747          ENDDO
3748
3749!
3750!--       Update the title attribute on file
3751!--       In order to avoid 'data mode' errors if updated attributes are larger
3752!--       than their original size, NF90_PUT_ATT is called in 'define mode'
3753!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
3754!--       performance loss due to data copying; an alternative strategy would be
3755!--       to ensure equal attribute size in a job chain. Maybe revise later.
3756          nc_stat = NF90_REDEF( id_set_ts )
3757          CALL handle_netcdf_error( 'netcdf', 439 )
3758          nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'title', &
3759                                  TRIM( run_description_header ) )
3760          CALL handle_netcdf_error( 'netcdf', 267 )
3761          nc_stat = NF90_ENDDEF( id_set_ts )
3762          CALL handle_netcdf_error( 'netcdf', 440 )
3763          message_string = 'NetCDF file for time series ' // & 
3764                           'from previous run found.' //     &
3765                           '&This file will be extended.'
3766          CALL message( 'define_netcdf_header', 'PA0259', 0, 0, 0, 6, 0 )
3767
3768
3769       CASE ( 'sp_new' )
3770
3771!
3772!--       Define some global attributes of the dataset
3773          IF ( averaging_interval_sp /= 0.0 )  THEN
3774             WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
3775                                                            averaging_interval_sp
3776             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title',  &
3777                                     TRIM( run_description_header ) // &
3778                                     TRIM( time_average_text ) )
3779             CALL handle_netcdf_error( 'netcdf', 268 )
3780
3781             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval_sp
3782             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'time_avg', &
3783                                     TRIM( time_average_text ) )
3784          ELSE
3785             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
3786                                     TRIM( run_description_header ) )
3787          ENDIF
3788          CALL handle_netcdf_error( 'netcdf', 269 )
3789
3790!
3791!--       Define time coordinate for spectra (unlimited dimension)
3792          nc_stat = NF90_DEF_DIM( id_set_sp, 'time', NF90_UNLIMITED, &
3793                                  id_dim_time_sp )
3794          CALL handle_netcdf_error( 'netcdf', 270 )
3795
3796          nc_stat = NF90_DEF_VAR( id_set_sp, 'time', NF90_DOUBLE, &
3797                                  id_dim_time_sp, id_var_time_sp )
3798          CALL handle_netcdf_error( 'netcdf', 271 )
3799
3800          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_time_sp, 'units', 'seconds')
3801          CALL handle_netcdf_error( 'netcdf', 272 )
3802
3803!
3804!--       Define the spatial dimensions and coordinates for spectra.
3805!--       First, determine the number of vertical levels for which spectra
3806!--       are to be output.
3807          ns = 1
3808          DO WHILE ( comp_spectra_level(ns) /= 999999  .AND.  ns <= 100 )
3809             ns = ns + 1
3810          ENDDO
3811          ns = ns - 1
3812
3813!
3814!--       Define vertical coordinate grid (zu grid)
3815          nc_stat = NF90_DEF_DIM( id_set_sp, 'zu_sp', ns, id_dim_zu_sp )
3816          CALL handle_netcdf_error( 'netcdf', 273 )
3817
3818          nc_stat = NF90_DEF_VAR( id_set_sp, 'zu_sp', NF90_DOUBLE, &
3819                                  id_dim_zu_sp, id_var_zu_sp )
3820          CALL handle_netcdf_error( 'netcdf', 274 )
3821
3822          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_zu_sp, 'units', 'meters' )
3823          CALL handle_netcdf_error( 'netcdf', 275 )
3824
3825!
3826!--       Define vertical coordinate grid (zw grid)
3827          nc_stat = NF90_DEF_DIM( id_set_sp, 'zw_sp', ns, id_dim_zw_sp )
3828          CALL handle_netcdf_error( 'netcdf', 276 )
3829
3830          nc_stat = NF90_DEF_VAR( id_set_sp, 'zw_sp', NF90_DOUBLE, &
3831                                  id_dim_zw_sp, id_var_zw_sp )
3832          CALL handle_netcdf_error( 'netcdf', 277 )
3833
3834          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_zw_sp, 'units', 'meters' )
3835          CALL handle_netcdf_error( 'netcdf', 278 )
3836
3837!
3838!--       Define x-axis
3839          nc_stat = NF90_DEF_DIM( id_set_sp, 'k_x', nx/2, id_dim_x_sp )
3840          CALL handle_netcdf_error( 'netcdf', 279 )
3841
3842          nc_stat = NF90_DEF_VAR( id_set_sp, 'k_x', NF90_DOUBLE, id_dim_x_sp, &
3843                                  id_var_x_sp )
3844          CALL handle_netcdf_error( 'netcdf', 280 )
3845
3846          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_x_sp, 'units', 'm-1' )
3847          CALL handle_netcdf_error( 'netcdf', 281 )
3848
3849!
3850!--       Define y-axis
3851          nc_stat = NF90_DEF_DIM( id_set_sp, 'k_y', ny/2, id_dim_y_sp )
3852          CALL handle_netcdf_error( 'netcdf', 282 )
3853
3854          nc_stat = NF90_DEF_VAR( id_set_sp, 'k_y', NF90_DOUBLE, id_dim_y_sp, &
3855                                  id_var_y_sp )
3856          CALL handle_netcdf_error( 'netcdf', 283 )
3857
3858          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_y_sp, 'units', 'm-1' )
3859          CALL handle_netcdf_error( 'netcdf', 284 )
3860
3861!
3862!--       Define the variables
3863          var_list = ';'
3864          i = 1
3865          DO WHILE ( data_output_sp(i) /= ' '  .AND.  i <= 10 )
3866!
3867!--          First check for the vertical grid
3868             found = .TRUE.
3869             SELECT CASE ( data_output_sp(i) )
3870!
3871!--             Most variables are defined on the zu levels
3872                CASE ( 'e', 'lpt', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c',  &
3873                       'ql_v', 'ql_vp', 'qv', 'rho', 's', 'sa', 'u', 'v', 'vpt' )
3874
3875                   grid_z = 'zu'
3876!
3877!--             zw levels
3878                CASE ( 'w' )
3879
3880                   grid_z = 'zw'
3881
3882                CASE DEFAULT
3883!
3884!--                Check for user-defined quantities (found, grid_x and grid_y
3885!--                are dummies)
3886                   CALL user_define_netcdf_grid( data_output_sp(i), found,  &
3887                                                 grid_x, grid_y, grid_z )
3888
3889             END SELECT
3890
3891             IF ( INDEX( spectra_direction(i), 'x' ) /= 0 )  THEN
3892
3893!
3894!--             Define the variable
3895                netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
3896                IF ( TRIM( grid_z ) == 'zw' )  THEN
3897                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
3898                                           nc_precision(7), (/ id_dim_x_sp, &
3899                                           id_dim_zw_sp, id_dim_time_sp /), &
3900                                           id_var_dospx(i) )
3901                ELSE
3902                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
3903                                           nc_precision(7), (/ id_dim_x_sp, &
3904                                           id_dim_zu_sp, id_dim_time_sp /), &
3905                                           id_var_dospx(i) )
3906                ENDIF
3907                CALL handle_netcdf_error( 'netcdf', 285 )
3908
3909                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospx(i), &
3910                                        'long_name', netcdf_var_name )
3911                CALL handle_netcdf_error( 'netcdf', 286 )
3912
3913                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospx(i), &
3914                                        'units', 'unknown' )
3915                CALL handle_netcdf_error( 'netcdf', 287 )
3916
3917                var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
3918
3919             ENDIF
3920
3921             IF ( INDEX( spectra_direction(i), 'y' ) /= 0 )  THEN
3922
3923!
3924!--             Define the variable
3925                netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
3926                IF ( TRIM( grid_z ) == 'zw' )  THEN
3927                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
3928                                           nc_precision(7), (/ id_dim_y_sp, &
3929                                           id_dim_zw_sp, id_dim_time_sp /), &
3930                                           id_var_dospy(i) )
3931                ELSE
3932                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
3933                                           nc_precision(7), (/ id_dim_y_sp, &
3934                                           id_dim_zu_sp, id_dim_time_sp /), &
3935                                           id_var_dospy(i) )
3936                ENDIF
3937                CALL handle_netcdf_error( 'netcdf', 288 )
3938
3939                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospy(i), &
3940                                        'long_name', netcdf_var_name )
3941                CALL handle_netcdf_error( 'netcdf', 289 )
3942
3943                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospy(i), &
3944                                        'units', 'unknown' )
3945                CALL handle_netcdf_error( 'netcdf', 290 )
3946
3947                var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
3948
3949             ENDIF
3950
3951             i = i + 1
3952
3953          ENDDO
3954
3955!
3956!--       Write the list of variables as global attribute (this is used by
3957!--       restart runs)
3958          nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', var_list )
3959          CALL handle_netcdf_error( 'netcdf', 291 )
3960
3961!
3962!--       Leave NetCDF define mode
3963          nc_stat = NF90_ENDDEF( id_set_sp )
3964          CALL handle_netcdf_error( 'netcdf', 292 )
3965
3966!
3967!--       Write axis data: zu_sp, zw_sp, k_x, k_y
3968          ALLOCATE( netcdf_data(1:ns) )
3969
3970!
3971!--       Write zu data
3972          netcdf_data(1:ns) = zu( comp_spectra_level(1:ns) )
3973          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zu_sp, netcdf_data, &
3974                                  start = (/ 1 /), count = (/ ns /) )
3975          CALL handle_netcdf_error( 'netcdf', 293 )
3976
3977!
3978!--       Write zw data
3979          netcdf_data(1:ns) = zw( comp_spectra_level(1:ns) )
3980          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zw_sp, netcdf_data, &
3981                                  start = (/ 1 /), count = (/ ns /) )
3982          CALL handle_netcdf_error( 'netcdf', 294 )
3983
3984          DEALLOCATE( netcdf_data )
3985
3986!
3987!--       Write data for x and y axis (wavenumbers)
3988          ALLOCATE( netcdf_data(nx/2) )
3989          DO  i = 1, nx/2
3990             netcdf_data(i) = 2.0 * pi * i / ( dx * ( nx + 1 ) )
3991          ENDDO
3992
3993          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_x_sp, netcdf_data, &
3994                                  start = (/ 1 /), count = (/ nx/2 /) )
3995          CALL handle_netcdf_error( 'netcdf', 295 )
3996
3997          DEALLOCATE( netcdf_data )
3998
3999          ALLOCATE( netcdf_data(ny/2) )
4000          DO  i = 1, ny/2
4001             netcdf_data(i) = 2.0 * pi * i / ( dy * ( ny + 1 ) )
4002          ENDDO
4003
4004          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_y_sp, netcdf_data, &
4005                                  start = (/ 1 /), count = (/ ny/2 /) )
4006          CALL handle_netcdf_error( 'netcdf', 296 )
4007
4008          DEALLOCATE( netcdf_data )
4009
4010
4011       CASE ( 'sp_ext' )
4012
4013!
4014!--       Get the list of variables and compare with the actual run.
4015!--       First var_list_old has to be reset, since GET_ATT does not assign
4016!--       trailing blanks.
4017          var_list_old = ' '
4018          nc_stat = NF90_GET_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', &
4019                                  var_list_old )
4020          CALL handle_netcdf_error( 'netcdf', 297 )
4021          var_list = ';'
4022          i = 1
4023          DO WHILE ( data_output_sp(i) /= ' '  .AND.  i <= 10 )
4024
4025             IF ( INDEX( spectra_direction(i), 'x' ) /= 0 )  THEN
4026                netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
4027                var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
4028             ENDIF
4029
4030             IF ( INDEX( spectra_direction(i), 'y' ) /= 0 )  THEN
4031                netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
4032                var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
4033             ENDIF
4034
4035             i = i + 1
4036
4037          ENDDO
4038
4039          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
4040             message_string = 'NetCDF file for spectra  ' //                 &
4041                              'from previous run found,' //                  &
4042                              '& but this file cannot be extended due to' // &
4043                              ' variable mismatch.' //                       &
4044                              '&New file is created instead.'
4045             CALL message( 'define_netcdf_header', 'PA0260', 0, 1, 0, 6, 0 )
4046             extend = .FALSE.
4047             RETURN
4048          ENDIF
4049
4050!
4051!--       Determine the number of current vertical levels for which spectra
4052!--       shall be output
4053          ns = 1
4054          DO WHILE ( comp_spectra_level(ns) /= 999999  .AND.  ns <= 100 )
4055             ns = ns + 1
4056          ENDDO
4057          ns = ns - 1
4058
4059!
4060!--       Get and compare the number of vertical levels
4061          nc_stat = NF90_INQ_VARID( id_set_sp, 'zu_sp', id_var_zu_sp )
4062          CALL handle_netcdf_error( 'netcdf', 298 )
4063
4064          nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_zu_sp, &
4065                                           dimids = id_dim_zu_sp_old )
4066          CALL handle_netcdf_error( 'netcdf', 299 )
4067          id_dim_zu_sp = id_dim_zu_sp_old(1)
4068
4069          nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_zu_sp, &
4070                                            len = ns_old )
4071          CALL handle_netcdf_error( 'netcdf', 300 )
4072
4073          IF ( ns /= ns_old )  THEN
4074             message_string = 'NetCDF file for spectra ' //                 &
4075                              ' from previous run found,' //                &
4076                              '&but this file cannot be extended due to' // &
4077                              ' mismatch in number of' //                   &
4078                              '&vertical levels.' //                        &
4079                              '&New file is created instead.'
4080             CALL message( 'define_netcdf_header', 'PA0261', 0, 1, 0, 6, 0 )
4081             extend = .FALSE.
4082             RETURN
4083          ENDIF
4084
4085!
4086!--       Get and compare the heights of the cross sections
4087          ALLOCATE( netcdf_data(1:ns_old) )
4088
4089          nc_stat = NF90_GET_VAR( id_set_sp, id_var_zu_sp, netcdf_data )
4090          CALL handle_netcdf_error( 'netcdf', 301 )
4091
4092          DO  i = 1, ns
4093             IF ( zu(comp_spectra_level(i)) /= netcdf_data(i) )  THEN
4094                message_string = 'NetCDF file for spectra ' //                 &
4095                                 ' from previous run found,' //                &
4096                                 '&but this file cannot be extended due to' // &
4097                                 ' mismatch in heights of' //                  &
4098                                 '&vertical levels.' //                        &
4099                                 '&New file is created instead.'
4100                CALL message( 'define_netcdf_header', 'PA0262', 0, 1, 0, 6, 0 )
4101                extend = .FALSE.
4102                RETURN
4103             ENDIF
4104          ENDDO
4105
4106          DEALLOCATE( netcdf_data )
4107
4108!
4109!--       Get the id of the time coordinate (unlimited coordinate) and its
4110!--       last index on the file. The next time level is plsp..count+1.
4111!--       The current time must be larger than the last output time
4112!--       on the file.
4113          nc_stat = NF90_INQ_VARID( id_set_sp, 'time', id_var_time_sp )
4114          CALL handle_netcdf_error( 'netcdf', 302 )
4115
4116          nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_time_sp, &
4117                                           dimids = id_dim_time_old )
4118          CALL handle_netcdf_error( 'netcdf', 303 )
4119          id_dim_time_sp = id_dim_time_old(1)
4120
4121          nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_time_sp, &
4122                                            len = dosp_time_count )
4123          CALL handle_netcdf_error( 'netcdf', 304 )
4124
4125          nc_stat = NF90_GET_VAR( id_set_sp, id_var_time_sp,        &
4126                                  last_time_coordinate,             &
4127                                  start = (/ dosp_time_count /), &
4128                                  count = (/ 1 /) )
4129          CALL handle_netcdf_error( 'netcdf', 305 )
4130
4131          IF ( last_time_coordinate(1) >= simulated_time )  THEN
4132             message_string = 'NetCDF file for spectra ' //                 &
4133                              'from previous run found,' //                 &
4134                              '&but this file cannot be extended becaus' // &
4135                              'e the current output time' //                & 
4136                              '&is less or equal than the last output t' // &
4137                              'ime on this file.' //                        &
4138                              '&New file is created instead.'
4139             CALL message( 'define_netcdf_header', 'PA0263', 0, 1, 0, 6, 0 )
4140             dosp_time_count = 0
4141             extend = .FALSE.
4142             RETURN
4143          ENDIF
4144
4145!
4146!--       Dataset seems to be extendable.
4147!--       Now get the variable ids.
4148          i = 1
4149          DO WHILE ( data_output_sp(i) /= ' '  .AND.  i <= 10 )
4150
4151             IF ( INDEX( spectra_direction(i), 'x' ) /= 0 )  THEN
4152                netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
4153                nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
4154                                          id_var_dospx(i) )
4155                CALL handle_netcdf_error( 'netcdf', 306 )
4156             ENDIF
4157
4158             IF ( INDEX( spectra_direction(i), 'y' ) /= 0 )  THEN
4159                netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
4160                nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
4161                                          id_var_dospy(i) )
4162                CALL handle_netcdf_error( 'netcdf', 307 )
4163             ENDIF
4164
4165             i = i + 1
4166
4167          ENDDO
4168
4169!
4170!--       Update the title attribute on file
4171!--       In order to avoid 'data mode' errors if updated attributes are larger
4172!--       than their original size, NF90_PUT_ATT is called in 'define mode'
4173!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
4174!--       performance loss due to data copying; an alternative strategy would be
4175!--       to ensure equal attribute size in a job chain. Maybe revise later.
4176          nc_stat = NF90_REDEF( id_set_sp )
4177          CALL handle_netcdf_error( 'netcdf', 441 )
4178          IF ( averaging_interval_sp /= 0.0 )  THEN
4179             WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
4180                                                           averaging_interval_sp
4181             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title',  &
4182                                     TRIM( run_description_header ) // &
4183                                     TRIM( time_average_text ) )
4184             CALL handle_netcdf_error( 'netcdf', 308 )
4185
4186             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval_sp
4187             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'time_avg', &
4188                                     TRIM( time_average_text ) )
4189          ELSE
4190             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
4191                                     TRIM( run_description_header ) )
4192          ENDIF
4193          CALL handle_netcdf_error( 'netcdf', 309 )
4194          nc_stat = NF90_ENDDEF( id_set_sp )
4195          CALL handle_netcdf_error( 'netcdf', 442 )
4196          message_string = 'NetCDF file for spectra ' //     & 
4197                           'from previous run found.' //     &
4198                           '&This file will be extended.'
4199          CALL message( 'define_netcdf_header', 'PA0264', 0, 0, 0, 6, 0 )
4200
4201
4202       CASE ( 'pt_new' )
4203
4204!
4205!--       Define some global attributes of the dataset
4206          nc_stat = NF90_PUT_ATT( id_set_prt, NF90_GLOBAL, 'title', &
4207                                  TRIM( run_description_header ) )
4208          CALL handle_netcdf_error( 'netcdf', 310 )
4209
4210!
4211!--       Define time coordinate for particles (unlimited dimension)
4212          nc_stat = NF90_DEF_DIM( id_set_prt, 'time', NF90_UNLIMITED, &
4213                                  id_dim_time_prt )
4214          CALL handle_netcdf_error( 'netcdf', 311 )
4215
4216          nc_stat = NF90_DEF_VAR( id_set_prt, 'time', NF90_DOUBLE, &
4217                                  id_dim_time_prt, id_var_time_prt )
4218          CALL handle_netcdf_error( 'netcdf', 312 )
4219
4220          nc_stat = NF90_PUT_ATT( id_set_prt, id_var_time_prt, 'units', &
4221                                  'seconds' )
4222          CALL handle_netcdf_error( 'netcdf', 313 )
4223
4224!
4225!--       Define particle coordinate (maximum particle number)
4226          IF ( netcdf_data_format < 3 )  THEN
4227             nc_stat = NF90_DEF_DIM( id_set_prt, 'particle_number', &
4228                                     maximum_number_of_particles, id_dim_prtnum)
4229          ELSE
4230!
4231!--          NetCDF4 allows more than one unlimited dimension
4232             nc_stat = NF90_DEF_DIM( id_set_prt, 'particle_number', &
4233                                     NF90_UNLIMITED, id_dim_prtnum)
4234          ENDIF
4235          CALL handle_netcdf_error( 'netcdf', 314 )
4236
4237          nc_stat = NF90_DEF_VAR( id_set_prt, 'particle_number', NF90_DOUBLE, &
4238                                  id_dim_prtnum, id_var_prtnum )
4239          CALL handle_netcdf_error( 'netcdf', 315 )
4240
4241          nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prtnum, 'units', &
4242                                  'particle number' )
4243          CALL handle_netcdf_error( 'netcdf', 316 )
4244
4245!
4246!--       Define variable which contains the real number of particles in use
4247          nc_stat = NF90_DEF_VAR( id_set_prt, 'real_num_of_prt', NF90_INT, &
4248                                  id_dim_time_prt, id_var_rnop_prt )
4249          CALL handle_netcdf_error( 'netcdf', 317 )
4250
4251          nc_stat = NF90_PUT_ATT( id_set_prt, id_var_rnop_prt, 'units', &
4252                                  'particle number' )
4253          CALL handle_netcdf_error( 'netcdf', 318 )
4254
4255!
4256!--       Define the variables
4257          DO  i = 1, 17
4258
4259             nc_stat = NF90_DEF_VAR( id_set_prt, prt_var_names(i),         &
4260                                     nc_precision(8),                      &
4261                                     (/ id_dim_prtnum, id_dim_time_prt /), &
4262                                     id_var_prt(i) )
4263             CALL handle_netcdf_error( 'netcdf', 319 )
4264
4265             nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prt(i), &
4266                                     'long_name', TRIM( prt_var_names(i) ) )
4267             CALL handle_netcdf_error( 'netcdf', 320 )
4268
4269             nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prt(i), &
4270                                     'units', TRIM( prt_var_units(i) ) )
4271             CALL handle_netcdf_error( 'netcdf', 321 )
4272
4273          ENDDO
4274
4275!
4276!--       Leave NetCDF define mode
4277          nc_stat = NF90_ENDDEF( id_set_prt )
4278          CALL handle_netcdf_error( 'netcdf', 322 )
4279
4280
4281       CASE ( 'pt_ext' )
4282
4283!
4284!--       Get the id of the time coordinate (unlimited coordinate) and its
4285!--       last index on the file. The next time level is prt..count+1.
4286!--       The current time must be larger than the last output time
4287!--       on the file.
4288          nc_stat = NF90_INQ_VARID( id_set_prt, 'time', id_var_time_prt )
4289          CALL handle_netcdf_error( 'netcdf', 323 )
4290
4291          nc_stat = NF90_INQUIRE_VARIABLE( id_set_prt, id_var_time_prt, &
4292                                           dimids = id_dim_time_old )
4293          CALL handle_netcdf_error( 'netcdf', 324 )
4294          id_dim_time_prt = id_dim_time_old(1)
4295
4296          nc_stat = NF90_INQUIRE_DIMENSION( id_set_prt, id_dim_time_prt, &
4297                                            len = prt_time_count )
4298          CALL handle_netcdf_error( 'netcdf', 325 )
4299
4300          nc_stat = NF90_GET_VAR( id_set_prt, id_var_time_prt,  &
4301                                  last_time_coordinate,         &
4302                                  start = (/ prt_time_count /), &
4303                                  count = (/ 1 /) )
4304          CALL handle_netcdf_error( 'netcdf', 326 )
4305
4306          IF ( last_time_coordinate(1) >= simulated_time )  THEN
4307             message_string = 'NetCDF file for particles ' //               &
4308                              'from previous run found,' //                 &
4309                              '&but this file cannot be extended becaus' // &
4310                              'e the current output time' //                &
4311                              '&is less or equal than the last output t' // &
4312                              'ime on this file.' //                        &
4313                              '&New file is created instead.'
4314             CALL message( 'define_netcdf_header', 'PA0265', 0, 1, 0, 6, 0 )
4315             prt_time_count = 0
4316             extend = .FALSE.
4317             RETURN
4318          ENDIF
4319
4320!
4321!--       Dataset seems to be extendable.
4322!--       Now get the variable ids.
4323          nc_stat = NF90_INQ_VARID( id_set_prt, 'real_num_of_prt', &
4324                                    id_var_rnop_prt )
4325          CALL handle_netcdf_error( 'netcdf', 327 )
4326
4327          DO  i = 1, 17
4328
4329             nc_stat = NF90_INQ_VARID( id_set_prt, prt_var_names(i), &
4330                                       id_var_prt(i) )
4331             CALL handle_netcdf_error( 'netcdf', 328 )
4332
4333          ENDDO
4334
4335          message_string = 'NetCDF file for particles ' // & 
4336                           'from previous run found.' //   &
4337                           '&This file will be extended.'
4338          CALL message( 'define_netcdf_header', 'PA0266', 0, 0, 0, 6, 0 )
4339       
4340
4341
4342       CASE ( 'ps_new' )
4343
4344!
4345!--       Define some global attributes of the dataset
4346          nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
4347                                  TRIM( run_description_header ) )
4348          CALL handle_netcdf_error( 'netcdf', 396 )
4349
4350!
4351!--       Define time coordinate for particle time series (unlimited dimension)
4352          nc_stat = NF90_DEF_DIM( id_set_pts, 'time', NF90_UNLIMITED, &
4353                                  id_dim_time_pts )
4354          CALL handle_netcdf_error( 'netcdf', 397 )
4355
4356          nc_stat = NF90_DEF_VAR( id_set_pts, 'time', NF90_DOUBLE, &
4357                                  id_dim_time_pts, id_var_time_pts )
4358          CALL handle_netcdf_error( 'netcdf', 398 )
4359
4360          nc_stat = NF90_PUT_ATT( id_set_pts, id_var_time_pts, 'units', &
4361                                  'seconds')
4362          CALL handle_netcdf_error( 'netcdf', 399 )
4363
4364!
4365!--       Define the variables. If more than one particle group is defined,
4366!--       define seperate variables for each group
4367          var_list = ';'
4368          DO  i = 1, dopts_num
4369
4370             DO  j = 0, number_of_particle_groups
4371
4372                IF ( j == 0 )  THEN
4373                   suffix1 = ''
4374                ELSE
4375                   WRITE ( suffix1, '(''_'',I2.2)' )  j
4376                ENDIF
4377
4378                nc_stat = NF90_DEF_VAR( id_set_pts,                         &
4379                                        TRIM( dopts_label(i) ) // suffix1,  &
4380                                        nc_precision(6),                    &
4381                                        (/ id_dim_time_pts /),              &
4382                                        id_var_dopts(i,j) )
4383                CALL handle_netcdf_error( 'netcdf', 400 )
4384
4385                IF ( j == 0 )  THEN
4386                   nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
4387                                           'long_name', &
4388                                           TRIM( dopts_label(i) ) )
4389                ELSE
4390                   nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
4391                                           'long_name', &
4392                                           TRIM( dopts_label(i) ) // ' PG ' // &
4393                                           suffix1(2:3) )
4394                ENDIF
4395                CALL handle_netcdf_error( 'netcdf', 401 )
4396
4397                nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
4398                                        'units', TRIM( dopts_unit(i) ) )
4399                CALL handle_netcdf_error( 'netcdf', 402 )
4400
4401                var_list = TRIM( var_list ) // TRIM( dopts_label(i) ) // &
4402                           suffix1 // ';'
4403
4404                IF ( number_of_particle_groups == 1 )  EXIT
4405
4406             ENDDO
4407
4408          ENDDO
4409
4410!
4411!--       Write the list of variables as global attribute (this is used by
4412!--       restart runs)
4413          nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
4414                                  var_list )
4415          CALL handle_netcdf_error( 'netcdf', 403 )
4416
4417
4418!
4419!--       Leave NetCDF define mode
4420          nc_stat = NF90_ENDDEF( id_set_pts )
4421          CALL handle_netcdf_error( 'netcdf', 404 )
4422
4423
4424       CASE ( 'ps_ext' )
4425
4426!
4427!--       Get the list of variables and compare with the actual run.
4428!--       First var_list_old has to be reset, since GET_ATT does not assign
4429!--       trailing blanks.
4430          var_list_old = ' '
4431          nc_stat = NF90_GET_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
4432                                  var_list_old )
4433          CALL handle_netcdf_error( 'netcdf', 405 )
4434
4435          var_list = ';'
4436          i = 1
4437          DO  i = 1, dopts_num
4438
4439             DO  j = 0, number_of_particle_groups
4440
4441                IF ( j == 0 )  THEN
4442                   suffix1 = ''
4443                ELSE
4444                   WRITE ( suffix1, '(''_'',I2.2)' )  j
4445                ENDIF
4446
4447                var_list = TRIM( var_list ) // TRIM( dopts_label(i) ) // &
4448                           suffix1 // ';'
4449
4450                IF ( number_of_particle_groups == 1 )  EXIT
4451
4452             ENDDO
4453
4454          ENDDO
4455
4456          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
4457             message_string = 'NetCDF file for particle time series ' //     &
4458                              'from previous run found,' //                  &
4459                              '& but this file cannot be extended due to' // &
4460                              ' variable mismatch.' //                       &
4461                              '&New file is created instead.'
4462             CALL message( 'define_netcdf_header', 'PA0267', 0, 1, 0, 6, 0 )
4463             extend = .FALSE.
4464             RETURN
4465          ENDIF
4466
4467!
4468!--       Get the id of the time coordinate (unlimited coordinate) and its
4469!--       last index on the file. The next time level is dots..count+1.
4470!--       The current time must be larger than the last output time
4471!--       on the file.
4472          nc_stat = NF90_INQ_VARID( id_set_pts, 'time', id_var_time_pts )
4473          CALL handle_netcdf_error( 'netcdf', 406 )
4474
4475          nc_stat = NF90_INQUIRE_VARIABLE( id_set_pts, id_var_time_pts, &
4476                                           dimids = id_dim_time_old )
4477          CALL handle_netcdf_error( 'netcdf', 407 )
4478          id_dim_time_pts = id_dim_time_old(1)
4479
4480          nc_stat = NF90_INQUIRE_DIMENSION( id_set_pts, id_dim_time_pts, &
4481                                            len = dopts_time_count )
4482          CALL handle_netcdf_error( 'netcdf', 408 )
4483
4484          nc_stat = NF90_GET_VAR( id_set_pts, id_var_time_pts,    &
4485                                  last_time_coordinate,           &
4486                                  start = (/ dopts_time_count /), &
4487                                  count = (/ 1 /) )
4488          CALL handle_netcdf_error( 'netcdf', 409 )
4489
4490          IF ( last_time_coordinate(1) >= simulated_time )  THEN
4491             message_string = 'NetCDF file for particle time series ' //    &
4492                              'from previous run found,' //                 &
4493                              '&but this file cannot be extended becaus' // &
4494                              'e the current output time' //                &
4495                              '&is less or equal than the last output t' // &
4496                              'ime on this file.' //                        &
4497                              '&New file is created instead.'
4498             CALL message( 'define_netcdf_header', 'PA0268', 0, 1, 0, 6, 0 )
4499             dopts_time_count = 0
4500             extend = .FALSE.
4501             RETURN
4502          ENDIF
4503
4504!
4505!--       Dataset seems to be extendable.
4506!--       Now get the variable ids
4507          i = 1
4508          DO  i = 1, dopts_num
4509 
4510             DO  j = 0, number_of_particle_groups
4511
4512                IF ( j == 0 )  THEN
4513                   suffix1 = ''
4514                ELSE
4515                   WRITE ( suffix1, '(''_'',I2.2)' )  j
4516                ENDIF
4517
4518                netcdf_var_name = TRIM( dopts_label(i) ) // suffix1
4519
4520                nc_stat = NF90_INQ_VARID( id_set_pts, netcdf_var_name, &
4521                                          id_var_dopts(i,j) )
4522                CALL handle_netcdf_error( 'netcdf', 410 )
4523
4524                IF ( number_of_particle_groups == 1 )  EXIT
4525
4526             ENDDO
4527
4528          ENDDO
4529
4530!
4531!--       Update the title attribute on file
4532!--       In order to avoid 'data mode' errors if updated attributes are larger
4533!--       than their original size, NF90_PUT_ATT is called in 'define mode'
4534!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
4535!--       performance loss due to data copying; an alternative strategy would be
4536!--       to ensure equal attribute size in a job chain. Maybe revise later.
4537          nc_stat = NF90_REDEF( id_set_pts )
4538          CALL handle_netcdf_error( 'netcdf', 443 )
4539          nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
4540                                  TRIM( run_description_header ) )
4541          CALL handle_netcdf_error( 'netcdf', 411 )
4542          nc_stat = NF90_ENDDEF( id_set_pts )
4543          CALL handle_netcdf_error( 'netcdf', 444 )
4544          message_string = 'NetCDF file for particle time series ' // & 
4545                           'from previous run found.' //              &
4546                           '&This file will be extended.'
4547          CALL message( 'define_netcdf_header', 'PA0269', 0, 0, 0, 6, 0 )
4548
4549
4550       CASE DEFAULT
4551
4552          message_string = 'mode "' // TRIM( mode) // '" not supported'
4553          CALL message( 'define_netcdf_header', 'PA0270', 0, 0, 0, 6, 0 )
4554
4555    END SELECT
4556
4557#endif
4558 END SUBROUTINE define_netcdf_header
4559
4560
4561SUBROUTINE handle_netcdf_error( routine_name, errno )
4562#if defined( __netcdf )
4563
4564!------------------------------------------------------------------------------!
4565!
4566! Description:
4567! ------------
4568! Prints out a text message corresponding to the current status.
4569!------------------------------------------------------------------------------!
4570
4571    USE control_parameters
4572    USE netcdf
4573    USE netcdf_control
4574    USE pegrid
4575
4576    IMPLICIT NONE
4577
4578    CHARACTER(LEN=6) ::  message_identifier
4579    CHARACTER(LEN=*) ::  routine_name
4580
4581    INTEGER ::  errno
4582
4583    IF ( nc_stat /= NF90_NOERR )  THEN
4584
4585       WRITE( message_identifier, '(''NC'',I4.4)' )  errno
4586       message_string = TRIM( NF90_STRERROR( nc_stat ) )
4587
4588       CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
4589
4590    ENDIF
4591
4592#endif
4593 END SUBROUTINE handle_netcdf_error
Note: See TracBrowser for help on using the repository browser.