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

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

two-moment cloud physics implemented

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