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

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

unused variables remove from several routines

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