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

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