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

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

last commit documented

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