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

Last change on this file since 1035 was 1035, checked in by raasch, 9 years ago

revisions r1031 and r1034 documented

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