source: palm/tags/release-3.2b/SOURCE/netcdf.f90 @ 4901

Last change on this file since 4901 was 77, checked in by raasch, 17 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:


General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:


Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

  • Property svn:keywords set to Id
File size: 142.5 KB
Line 
1#if defined( __ibmy_special )
2@PROCESS NOOPTimize
3#endif
4 SUBROUTINE define_netcdf_header( callmode, extend, av )
5
6!------------------------------------------------------------------------------!
7! Current revisions:
8! ------------------
9!
10!
11! Former revisions:
12! -----------------
13! $Id: netcdf.f90 77 2007-03-29 04:26:56Z banzhafs $
14!
15! 48 2007-03-06 12:28:36Z raasch
16! Output topography height information (zu_s_inner, zw_s_inner) to 2d-xy and 3d
17! datasets
18!
19! RCS Log replace by Id keyword, revision history cleaned up
20!
21! Revision 1.12  2006/09/26 19:35:16  raasch
22! Bugfix yv coordinates for yz cross sections
23!
24! Revision 1.1  2005/05/18 15:37:16  raasch
25! Initial revision
26!
27!
28! Description:
29! ------------
30! In case of extend = .FALSE.:
31! Define all necessary dimensions, axes and variables for the different
32! NetCDF datasets. This subroutine is called from check_open after a new
33! dataset is created. It leaves the open NetCDF files ready to write.
34!
35! In case of extend = .TRUE.:
36! Find out if dimensions and variables of an existing file match the values
37! of the actual run. If so, get all necessary informations (ids, etc.) from
38! this file.
39!
40! Parameter av can assume values 0 (non-averaged data) and 1 (time averaged
41! data)
42!------------------------------------------------------------------------------!
43#if defined( __netcdf )
44
45    USE arrays_3d
46    USE constants
47    USE control_parameters
48    USE grid_variables
49    USE indices
50    USE netcdf_control
51    USE pegrid
52    USE particle_attributes
53    USE profil_parameter
54    USE spectrum
55    USE statistics
56
57
58    IMPLICIT NONE
59
60    CHARACTER (LEN=2)              ::  suffix
61    CHARACTER (LEN=2), INTENT (IN) ::  callmode
62    CHARACTER (LEN=3)              ::  suffix1
63    CHARACTER (LEN=4)              ::  grid_x, grid_y, grid_z
64    CHARACTER (LEN=6)              ::  mode
65    CHARACTER (LEN=10)             ::  netcdf_var_name, netcdf_var_name_base, &
66                                       precision, var
67    CHARACTER (LEN=80)             ::  time_average_text
68    CHARACTER (LEN=2000)           ::  var_list, var_list_old
69
70    INTEGER ::  av, i, id_x, id_y, id_z, j, ns, ns_old, nz_old
71
72    INTEGER, DIMENSION(1) ::  id_dim_time_old, id_dim_x_yz_old,  &
73                              id_dim_y_xz_old, id_dim_zu_sp_old, &
74                              id_dim_zu_xy_old, id_dim_zu_3d_old
75
76    LOGICAL ::  found
77
78    LOGICAL, INTENT (INOUT) ::  extend
79
80    LOGICAL, SAVE ::  init_netcdf = .FALSE.
81
82    REAL, DIMENSION(1) ::  last_time_coordinate
83
84    REAL, DIMENSION(:), ALLOCATABLE ::  netcdf_data
85
86
87!
88!-- Initializing actions (return to calling routine check_parameters afterwards)
89    IF ( .NOT. init_netcdf )  THEN
90!
91!--    Check and set accuracy for NetCDF output. First set default value
92       nc_precision = NF90_REAL4
93
94       i = 1
95       DO  WHILE ( netcdf_precision(i) /= ' ' )
96          j = INDEX( netcdf_precision(i), '_' )
97          IF ( j == 0 )  THEN
98             IF ( myid == 0 )  THEN
99                PRINT*, '+++ define_netcdf_header: netcdf_precision must ',  &
100                             'contain a "_"    netcdf_precision(', i, ')="', &
101                             TRIM( netcdf_precision(i) ),'"'
102             ENDIF
103             CALL local_stop
104          ENDIF
105
106          var       = netcdf_precision(i)(1:j-1)
107          precision = netcdf_precision(i)(j+1:)
108
109          IF ( precision == 'NF90_REAL4' )  THEN
110             j = NF90_REAL4
111          ELSEIF ( precision == 'NF90_REAL8' )  THEN
112             j = NF90_REAL8
113          ELSE
114             IF ( myid == 0 )  THEN
115                PRINT*, '+++ define_netcdf_header: illegal netcdf precision: ',&
116                             'netcdf_precision(', i, ')="', &
117                             TRIM( netcdf_precision(i) ),'"'
118             ENDIF
119             CALL local_stop
120          ENDIF
121
122          SELECT CASE ( var )
123             CASE ( 'xy' )
124                nc_precision(1) = j
125             CASE ( 'xz' )
126                nc_precision(2) = j
127             CASE ( 'yz' )
128                nc_precision(3) = j
129             CASE ( '2d' )
130                nc_precision(1:3) = j
131             CASE ( '3d' )
132                nc_precision(4) = j
133             CASE ( 'pr' )
134                nc_precision(5) = j
135             CASE ( 'ts' )
136                nc_precision(6) = j
137             CASE ( 'sp' )
138                nc_precision(7) = j
139             CASE ( 'prt' )
140                nc_precision(8) = j
141             CASE ( 'all' )
142                nc_precision    = j
143
144             CASE DEFAULT
145                IF ( myid == 0 )  THEN
146                   PRINT*, '+++ define_netcdf_header: unknown variable in ',   &
147                              'inipar assignment: netcdf_precision(', i, ')="',&
148                                TRIM( netcdf_precision(i) ),'"'
149                ENDIF
150                CALL local_stop               
151
152          END SELECT
153
154          i = i + 1
155          IF ( i > 10 )  EXIT
156       ENDDO
157
158       init_netcdf = .TRUE.
159
160       RETURN
161
162    ENDIF
163
164!
165!-- Determine the mode to be processed
166    IF ( extend )  THEN
167       mode = callmode // '_ext'
168    ELSE
169       mode = callmode // '_new'
170    ENDIF
171
172!
173!-- Select the mode to be processed. Possibilities are xy, xz, yz, pr and ts.
174    SELECT CASE ( mode )
175
176       CASE ( '3d_new' )
177
178!
179!--       Define some global attributes of the dataset
180          nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'Conventions', &
181                                  'COARDS' )
182          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 62 )
183
184          IF ( av == 0 )  THEN
185             time_average_text = ' '
186          ELSE
187             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
188                                                            averaging_interval
189          ENDIF
190          nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'title', &
191                                  TRIM( run_description_header ) //    &
192                                  TRIM( time_average_text ) )
193          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 63 )
194          IF ( av == 1 )  THEN
195             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
196             nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'time_avg', &
197                                     TRIM( time_average_text ) )
198             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 63 )
199          ENDIF
200
201!
202!--       Define time coordinate for volume data (unlimited dimension)
203          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'time', NF90_UNLIMITED, &
204                                  id_dim_time_3d(av) )
205          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 64 )
206
207          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'time', NF90_DOUBLE, &
208                                  id_dim_time_3d(av), id_var_time_3d(av) )
209          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 65 )
210
211          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_time_3d(av), 'units', &
212                                  'seconds')
213          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 66 )
214
215!
216!--       Define spatial dimensions and coordinates:
217!--       Define vertical coordinate grid (zu grid)
218          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'zu_3d', nz_do3d-nzb+1, &
219                                  id_dim_zu_3d(av) )
220          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 67 )
221
222          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zu_3d', NF90_DOUBLE, &
223                                  id_dim_zu_3d(av), id_var_zu_3d(av) )
224          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 68 )
225
226          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zu_3d(av), 'units', &
227                                  'meters' )
228          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 69 )
229
230!
231!--       Define vertical coordinate grid (zw grid)
232          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'zw_3d', nz_do3d-nzb+1, &
233                                  id_dim_zw_3d(av) )
234          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 70 )
235
236          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zw_3d', NF90_DOUBLE, &
237                                  id_dim_zw_3d(av), id_var_zw_3d(av) )
238          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 71 )
239
240          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zw_3d(av), 'units', &
241                                  'meters' )
242          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 72 )
243
244!
245!--       Define x-axis (for scalar position)
246          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'x', nx+2, id_dim_x_3d(av) )
247          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 73 )
248
249          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'x', NF90_DOUBLE, &
250                                  id_dim_x_3d(av), id_var_x_3d(av) )
251          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 74 )
252
253          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_x_3d(av), 'units', &
254                                  'meters' )
255          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 75 )
256
257!
258!--       Define x-axis (for u position)
259          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'xu', nx+2, id_dim_xu_3d(av) )
260          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 358 )
261
262          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'xu', NF90_DOUBLE, &
263                                  id_dim_xu_3d(av), id_var_xu_3d(av) )
264          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 359 )
265
266          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_xu_3d(av), 'units', &
267                                  'meters' )
268          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 360 )
269
270!
271!--       Define y-axis (for scalar position)
272          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'y', ny+2, id_dim_y_3d(av) )
273          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 76 )
274
275          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'y', NF90_DOUBLE, &
276                                  id_dim_y_3d(av), id_var_y_3d(av) )
277          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 77 )
278
279          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_y_3d(av), 'units', &
280                                  'meters' )
281          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 78 )
282
283!
284!--       Define y-axis (for v position)
285          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'yv', ny+2, id_dim_yv_3d(av) )
286          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 361 )
287
288          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'yv', NF90_DOUBLE, &
289                                  id_dim_yv_3d(av), id_var_yv_3d(av) )
290          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 362 )
291
292          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_yv_3d(av), 'units', &
293                                  'meters' )
294          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 363 )
295
296!
297!--       In case of non-flat topography define 2d-arrays containing the height
298!--       informations
299          IF ( TRIM( topography ) /= 'flat' )  THEN
300!
301!--          Define zusi = zu(nzb_s_inner)
302             nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zusi', NF90_DOUBLE,     &
303                                     (/ id_dim_x_3d(av), id_dim_y_3d(av) /), &
304                                     id_var_zusi_3d(av) )
305             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 413 )
306             
307             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zusi_3d(av), &
308                                     'units', 'meters' )
309             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 414 )
310             
311             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zusi_3d(av), &
312                                     'long_name', 'zu(nzb_s_inner)' )
313             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 415 )
314
315!             
316!--          Define zwwi = zw(nzb_w_inner)
317             nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zwwi', NF90_DOUBLE,     &
318                                     (/ id_dim_x_3d(av), id_dim_y_3d(av) /), &
319                                     id_var_zwwi_3d(av) )
320             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 416 )
321             
322             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zwwi_3d(av), &
323                                     'units', 'meters' )
324             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 417 )
325             
326             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zwwi_3d(av), &
327                                     'long_name', 'zw(nzb_w_inner)' )
328             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 418 )
329
330          ENDIF             
331
332
333!
334!--       Define the variables
335          var_list = ';'
336          i = 1
337
338          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
339
340!
341!--          Check for the grid
342             found = .TRUE.
343             SELECT CASE ( do3d(av,i) )
344!
345!--             Most variables are defined on the scalar grid
346                CASE ( 'e', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c', 'ql_v', &
347                       'ql_vp', 'qv', 's', 'vpt' )
348
349                   grid_x = 'x'
350                   grid_y = 'y'
351                   grid_z = 'zu'
352!
353!--             u grid
354                CASE ( 'u' )
355
356                   grid_x = 'xu'
357                   grid_y = 'y'
358                   grid_z = 'zu'
359!
360!--             v grid
361                CASE ( 'v' )
362
363                   grid_x = 'x'
364                   grid_y = 'yv'
365                   grid_z = 'zu'
366!
367!--             w grid
368                CASE ( 'w' )
369
370                   grid_x = 'x'
371                   grid_y = 'y'
372                   grid_z = 'zw'
373
374                CASE DEFAULT
375!
376!--                Check for user-defined quantities
377                   CALL user_define_netcdf_grid( do3d(av,i), found, grid_x, &
378                                                 grid_y, grid_z )
379
380             END SELECT
381
382!
383!--          Select the respective dimension ids
384             IF ( grid_x == 'x' )  THEN
385                id_x = id_dim_x_3d(av)
386             ELSEIF ( grid_x == 'xu' )  THEN
387                id_x = id_dim_xu_3d(av)
388             ENDIF
389
390             IF ( grid_y == 'y' )  THEN
391                id_y = id_dim_y_3d(av)
392             ELSEIF ( grid_y == 'yv' )  THEN
393                id_y = id_dim_yv_3d(av)
394             ENDIF
395
396             IF ( grid_z == 'zu' )  THEN
397                id_z = id_dim_zu_3d(av)
398             ELSEIF ( grid_z == 'zw' )  THEN
399                id_z = id_dim_zw_3d(av)
400             ENDIF
401
402!
403!--          Define the grid
404             nc_stat = NF90_DEF_VAR( id_set_3d(av), do3d(av,i),                &
405                                     nc_precision(4),                          &
406                                   (/ id_x, id_y, id_z, id_dim_time_3d(av) /), &
407                                     id_var_do3d(av,i) )
408
409             IF ( .NOT. found )  THEN
410                PRINT*, '+++ define_netcdf_header: no grid defined for', &
411                             ' variable ', do3d(av,i)
412             ENDIF
413
414             var_list = TRIM( var_list ) // TRIM( do3d(av,i) ) // ';'
415
416             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 79 )
417!
418!--          Store the 'real' name of the variable (with *, for example)
419!--          in the long_name attribute. This is evaluated by Ferret,
420!--          for example.
421             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_do3d(av,i), &
422                                     'long_name', do3d(av,i) )
423             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 80 )
424!
425!--          Define the variable's unit
426             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_do3d(av,i), &
427                                     'units', TRIM( do3d_unit(av,i) ) )
428             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 357 )
429
430             i = i + 1
431
432          ENDDO
433
434!
435!--       No arrays to output
436          IF ( i == 1 )  RETURN
437
438!
439!--       Write the list of variables as global attribute (this is used by
440!--       restart runs and by combine_plot_fields)
441          nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
442                                  var_list )
443          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 81 )
444
445!
446!--       Leave NetCDF define mode
447          nc_stat = NF90_ENDDEF( id_set_3d(av) )
448          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 82 )
449
450!
451!--       Write data for x and xu axis (shifted by -dx/2)
452          ALLOCATE( netcdf_data(0:nx+1) )
453
454          DO  i = 0, nx+1
455             netcdf_data(i) = i * dx
456          ENDDO
457
458          nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_x_3d(av), netcdf_data, &
459                                  start = (/ 1 /), count = (/ nx+2 /) )
460          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 83 )
461
462          DO  i = 0, nx+1
463             netcdf_data(i) = ( i - 0.5 ) * dx
464          ENDDO
465
466          nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_xu_3d(av), &
467                                  netcdf_data, start = (/ 1 /),    &
468                                  count = (/ nx+2 /) )
469          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 385 )
470
471          DEALLOCATE( netcdf_data )
472
473!
474!--       Write data for y and yv axis (shifted by -dy/2)
475          ALLOCATE( netcdf_data(0:ny+1) )
476
477          DO  i = 0, ny+1
478             netcdf_data(i) = i * dy
479          ENDDO
480
481          nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_y_3d(av), netcdf_data, &
482                                  start = (/ 1 /), count = (/ ny+2 /))
483          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 84 )
484
485          DO  i = 0, ny+1
486             netcdf_data(i) = ( i - 0.5 ) * dy
487          ENDDO
488
489          nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_yv_3d(av), &
490                                  netcdf_data, start = (/ 1 /),    &
491                                  count = (/ ny+2 /))
492          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 387 )
493
494          DEALLOCATE( netcdf_data )
495
496!
497!--       Write zu and zw data (vertical axes)
498          nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zu_3d(av),    &
499                                  zu(nzb:nz_do3d), start = (/ 1 /), &
500                                  count = (/ nz_do3d-nzb+1 /) )
501          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 85 )
502
503          nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zw_3d(av),    &
504                                  zw(nzb:nz_do3d), start = (/ 1 /), &
505                                  count = (/ nz_do3d-nzb+1 /) )
506          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 86 )
507
508
509!
510!--       In case of non-flat topography write height information
511          IF ( TRIM( topography ) /= 'flat' )  THEN
512
513             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av), &
514                                     zu_s_inner(0:nx+1,0:ny+1), &
515                                     start = (/ 1, 1 /), &
516                                     count = (/ nx+2, ny+2 /) )
517             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 419 )
518
519             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av), &
520                                     zw_w_inner(0:nx+1,0:ny+1), &
521                                     start = (/ 1, 1 /), &
522                                     count = (/ nx+2, ny+2 /) )
523             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 420 )
524
525          ENDIF
526
527       CASE ( '3d_ext' )
528
529!
530!--       Get the list of variables and compare with the actual run.
531!--       First var_list_old has to be reset, since GET_ATT does not assign
532!--       trailing blanks.
533          var_list_old = ' '
534          nc_stat = NF90_GET_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
535                                  var_list_old )
536          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 87 )
537
538          var_list = ';'
539          i = 1
540          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
541             var_list = TRIM(var_list) // TRIM( do3d(av,i) ) // ';'
542             i = i + 1
543          ENDDO
544
545          IF ( av == 0 )  THEN
546             var = '(3d)'
547          ELSE
548             var = '(3d_av)'
549          ENDIF
550
551          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
552             PRINT*, '+++ WARNING: NetCDF file for volume data ' // &
553                                   TRIM( var ) // ' from previuos run found,'
554             PRINT*, '             but this file cannot be extended due to' // &
555                                   ' variable mismatch.'
556             PRINT*, '             New file is created instead.'
557             extend = .FALSE.
558             RETURN
559          ENDIF
560
561!
562!--       Get and compare the number of vertical gridpoints
563          nc_stat = NF90_INQ_VARID( id_set_3d(av), 'zu_3d', id_var_zu_3d(av) )
564          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 88 )
565
566          nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_zu_3d(av), &
567                                           dimids = id_dim_zu_3d_old )
568          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 89 )
569          id_dim_zu_3d(av) = id_dim_zu_3d_old(1)
570
571          nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_zu_3d(av), &
572                                            len = nz_old )
573          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 90 )
574
575          IF ( nz_do3d-nzb+1 /= nz_old )  THEN
576             PRINT*, '+++ WARNING: NetCDF file for volume data ' // &
577                                   TRIM( var ) // ' from previuos run found,'
578             PRINT*, '             but this file cannot be extended due to' // &
579                                   ' mismatch in number of'
580             PRINT*, '             vertical grid points (nz_do3d).'
581             PRINT*, '             New file is created instead.'
582             extend = .FALSE.
583             RETURN
584          ENDIF
585
586!
587!--       Get the id of the time coordinate (unlimited coordinate) and its
588!--       last index on the file. The next time level is pl3d..count+1.
589!--       The current time must be larger than the last output time
590!--       on the file.
591          nc_stat = NF90_INQ_VARID( id_set_3d(av), 'time', id_var_time_3d(av) )
592          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 91 )
593
594          nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_time_3d(av), &
595                                           dimids = id_dim_time_old )
596          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 92 )
597          id_dim_time_3d(av) = id_dim_time_old(1)
598
599          nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_time_3d(av), &
600                                            len = do3d_time_count(av) )
601          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 93 )
602
603          nc_stat = NF90_GET_VAR( id_set_3d(av), id_var_time_3d(av), &
604                                  last_time_coordinate,              &
605                                  start = (/ do3d_time_count(av) /), &
606                                  count = (/ 1 /) )
607          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 94 )
608
609          IF ( last_time_coordinate(1) >= simulated_time )  THEN
610             PRINT*, '+++ WARNING: NetCDF file for volume data ' // &
611                                   TRIM( var ) // ' from previuos run found,'
612             PRINT*, '             but this file cannot be extended becaus' // &
613                                   'e the current output time'
614             PRINT*, '             is less or equal than the last output t' // &
615                                   'ime on this file.'
616             PRINT*, '             New file is created instead.'
617             do3d_time_count(av) = 0
618             extend = .FALSE.
619             RETURN
620          ENDIF
621
622!
623!--       Dataset seems to be extendable.
624!--       Now get the variable ids.
625          i = 1
626          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
627             nc_stat = NF90_INQ_VARID( id_set_3d(av), TRIM( do3d(av,i) ), &
628                                       id_var_do3d(av,i) )
629             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 95 )
630             i = i + 1
631          ENDDO
632
633!
634!--       Change the titel attribute on file
635          nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'title', &
636                                  TRIM( run_description_header ) )
637          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 96 )
638          PRINT*, '*** NetCDF file for volume data ' // TRIM( var ) // &
639                       ' from previous run found.'
640          PRINT*, '    This file will be extended.'
641
642
643       CASE ( 'xy_new' )
644
645!
646!--       Define some global attributes of the dataset
647          nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'Conventions', &
648                                  'COARDS' )
649          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 97 )
650
651          IF ( av == 0 )  THEN
652             time_average_text = ' '
653          ELSE
654             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
655                                                            averaging_interval
656          ENDIF
657          nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'title', &
658                                  TRIM( run_description_header ) //    &
659                                  TRIM( time_average_text ) )
660          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 98 )
661          IF ( av == 1 )  THEN
662             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
663             nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'time_avg', &
664                                     TRIM( time_average_text ) )
665             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 98 )
666          ENDIF
667
668!
669!--       Define time coordinate for xy sections (unlimited dimension)
670          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'time', NF90_UNLIMITED, &
671                                  id_dim_time_xy(av) )
672          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 99 )
673
674          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'time', NF90_DOUBLE, &
675                                  id_dim_time_xy(av), id_var_time_xy(av) )
676          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 100 )
677
678          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_time_xy(av), 'units', &
679                                  'seconds')
680          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 101 )
681
682!
683!--       Define the spatial dimensions and coordinates for xy-sections.
684!--       First, determine the number of horizontal sections to be written.
685          IF ( section(1,1) == -9999 )  THEN
686             RETURN
687          ELSE
688             ns = 1
689             DO WHILE ( section(ns,1) /= -9999  .AND.  ns <= 100 )
690                ns = ns + 1
691             ENDDO
692             ns = ns - 1
693          ENDIF
694
695!
696!--       Define vertical coordinate grid (zu grid)
697          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'zu_xy', ns, id_dim_zu_xy(av) )
698          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 102 )
699
700          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zu_xy', NF90_DOUBLE, &
701                                  id_dim_zu_xy(av), id_var_zu_xy(av) )
702          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 103 )
703
704          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zu_xy(av), 'units', &
705                                  'meters' )
706          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 104 )
707
708!
709!--       Define vertical coordinate grid (zw grid)
710          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'zw_xy', ns, id_dim_zw_xy(av) )
711          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 105 )
712
713          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zw_xy', NF90_DOUBLE, &
714                                  id_dim_zw_xy(av), id_var_zw_xy(av) )
715          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 106 )
716
717          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zw_xy(av), 'units', &
718                                  'meters' )
719          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 107 )
720
721!
722!--       Define a pseudo vertical coordinate grid for the surface variables
723!--       u* and t* to store their height level
724          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'zu1_xy', 1, &
725                                  id_dim_zu1_xy(av) )
726          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 108 )
727
728          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zu1_xy', NF90_DOUBLE, &
729                                  id_dim_zu1_xy(av), id_var_zu1_xy(av) )
730          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 109 )
731
732          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zu1_xy(av), 'units', &
733                                  'meters' )
734          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 110 )
735
736!
737!--       Define a variable to store the layer indices of the horizontal cross
738!--       sections, too
739          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'ind_z_xy', NF90_DOUBLE, &
740                                  id_dim_zu_xy(av), id_var_ind_z_xy(av) )
741          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 111 )
742
743          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_ind_z_xy(av), 'units', &
744                                  'gridpoints')
745          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 112 )
746
747!
748!--       Define x-axis (for scalar position)
749          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'x', nx+2, id_dim_x_xy(av) )
750          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 113 )
751
752          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'x', NF90_DOUBLE, &
753                                  id_dim_x_xy(av), id_var_x_xy(av) )
754          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 114 )
755
756          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_x_xy(av), 'units', &
757                                  'meters' )
758          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 115 )
759
760!
761!--       Define x-axis (for u position)
762          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'xu', nx+2, id_dim_xu_xy(av) )
763          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 388 )
764
765          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'xu', NF90_DOUBLE, &
766                                  id_dim_xu_xy(av), id_var_xu_xy(av) )
767          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 389 )
768
769          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_xu_xy(av), 'units', &
770                                  'meters' )
771          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 390 )
772
773!
774!--       Define y-axis (for scalar position)
775          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'y', ny+2, id_dim_y_xy(av) )
776          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 116 )
777
778          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'y', NF90_DOUBLE, &
779                                  id_dim_y_xy(av), id_var_y_xy(av) )
780          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 117 )
781
782          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_y_xy(av), 'units', &
783                                  'meters' )
784          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 118 )
785
786!
787!--       Define y-axis (for scalar position)
788          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'yv', ny+2, id_dim_yv_xy(av) )
789          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 364 )
790
791          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'yv', NF90_DOUBLE, &
792                                  id_dim_yv_xy(av), id_var_yv_xy(av) )
793          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 365 )
794
795          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_yv_xy(av), 'units', &
796                                  'meters' )
797          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 366 )
798
799!
800!--       In case of non-flat topography define 2d-arrays containing the height
801!--       informations
802          IF ( TRIM( topography ) /= 'flat' )  THEN
803!
804!--          Define zusi = zu(nzb_s_inner)
805             nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zusi', NF90_DOUBLE, &
806                                     (/ id_dim_x_xy(av), id_dim_y_xy(av) /), &
807                                     id_var_zusi_xy(av) )
808             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 421 )
809             
810             nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zusi_xy(av), &
811                                     'units', 'meters' )
812             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 422 )
813             
814             nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zusi_xy(av), &
815                                     'long_name', 'zu(nzb_s_inner)' )
816             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 423 )
817
818!             
819!--          Define zwwi = zw(nzb_w_inner)
820             nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zwwi', NF90_DOUBLE, &
821                                     (/ id_dim_x_xy(av), id_dim_y_xy(av) /), &
822                                     id_var_zwwi_xy(av) )
823             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 424 )
824             
825             nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zwwi_xy(av), &
826                                     'units', 'meters' )
827             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 425 )
828             
829             nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zwwi_xy(av), &
830                                     'long_name', 'zw(nzb_w_inner)' )
831             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 426 )
832
833          ENDIF
834
835
836!
837!--       Define the variables
838          var_list = ';'
839          i = 1
840
841          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
842
843             IF ( INDEX( do2d(av,i), 'xy' ) /= 0 )  THEN
844!
845!--             If there is a star in the variable name (u* or t*), it is a
846!--             surface variable. Define it with id_dim_zu1_xy.
847                IF ( INDEX( do2d(av,i), '*' ) /= 0 )  THEN
848!
849!--                First, remove those characters not allowed by NetCDF
850                   netcdf_var_name = do2d(av,i)
851                   CALL clean_netcdf_varname( netcdf_var_name )
852
853                   nc_stat = NF90_DEF_VAR( id_set_xy(av), netcdf_var_name,     &
854                                           nc_precision(1),                    &
855                                           (/ id_dim_x_xy(av), id_dim_y_xy(av),&
856                                              id_dim_zu1_xy(av),               &
857                                              id_dim_time_xy(av) /),           &
858                                           id_var_do2d(av,i) )
859
860                   var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
861
862                ELSE
863
864!
865!--                Check for the grid
866                   found = .TRUE.
867                   SELECT CASE ( do2d(av,i) )
868!
869!--                   Most variables are defined on the zu grid
870                      CASE ( 'e_xy', 'p_xy', 'pc_xy', 'pr_xy', 'pt_xy', 'q_xy',&
871                             'ql_xy', 'ql_c_xy', 'ql_v_xy', 'ql_vp_xy',        &
872                             'qv_xy', 's_xy', 'vpt_xy' )
873
874                         grid_x = 'x'
875                         grid_y = 'y'
876                         grid_z = 'zu'
877!
878!--                   u grid
879                      CASE ( 'u_xy' )
880
881                         grid_x = 'xu'
882                         grid_y = 'y'
883                         grid_z = 'zu'
884!
885!--                   v grid
886                      CASE ( 'v_xy' )
887
888                         grid_x = 'x'
889                         grid_y = 'yv'
890                         grid_z = 'zu'
891!
892!--                   w grid
893                      CASE ( 'w_xy' )
894
895                         grid_x = 'x'
896                         grid_y = 'y'
897                         grid_z = 'zw'
898
899                      CASE DEFAULT
900!
901!--                      Check for user-defined quantities
902                         CALL user_define_netcdf_grid( do2d(av,i), found, &
903                                                       grid_x, grid_y, grid_z )
904
905                   END SELECT
906
907!
908!--                Select the respective dimension ids
909                   IF ( grid_x == 'x' )  THEN
910                      id_x = id_dim_x_xy(av)
911                   ELSEIF ( grid_x == 'xu' )  THEN
912                      id_x = id_dim_xu_xy(av)
913                   ENDIF
914
915                   IF ( grid_y == 'y' )  THEN
916                      id_y = id_dim_y_xy(av)
917                   ELSEIF ( grid_y == 'yv' )  THEN
918                      id_y = id_dim_yv_xy(av)
919                   ENDIF
920
921                   IF ( grid_z == 'zu' )  THEN
922                      id_z = id_dim_zu_xy(av)
923                   ELSEIF ( grid_z == 'zw' )  THEN
924                      id_z = id_dim_zw_xy(av)
925                   ENDIF
926
927!
928!--                Define the grid
929                   nc_stat = NF90_DEF_VAR( id_set_xy(av), do2d(av,i),          &
930                                           nc_precision(1),                    &
931                                   (/ id_x, id_y, id_z, id_dim_time_xy(av) /), &
932                                           id_var_do2d(av,i) )
933
934                   IF ( .NOT. found )  THEN
935                      PRINT*, '+++ define_netcdf_header: no grid defined ', &
936                                   'for variable ', do2d(av,i)
937                   ENDIF
938
939                   var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
940
941                ENDIF
942
943                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 119 )
944!
945!--             Store the 'real' name of the variable (with *, for example)
946!--             in the long_name attribute. This is evaluated by Ferret,
947!--             for example.
948                nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_do2d(av,i), &
949                                        'long_name', do2d(av,i) )
950                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 120 )
951!
952!--             Define the variable's unit
953                nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_do2d(av,i), &
954                                        'units', TRIM( do2d_unit(av,i) ) )
955                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 354 )
956             ENDIF
957
958             i = i + 1
959
960          ENDDO
961
962!
963!--       No arrays to output. Close the netcdf file and return.
964          IF ( i == 1 )  RETURN
965
966!
967!--       Write the list of variables as global attribute (this is used by
968!--       restart runs and by combine_plot_fields)
969          nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'VAR_LIST', &
970                                  var_list )
971          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 121 )
972
973!
974!--       Leave NetCDF define mode
975          nc_stat = NF90_ENDDEF( id_set_xy(av) )
976          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 122 )
977
978!
979!--       Write axis data: z_xy, x, y
980          ALLOCATE( netcdf_data(1:ns) )
981
982!
983!--       Write zu data
984          DO  i = 1, ns
985             IF( section(i,1) == -1 )  THEN
986                netcdf_data(i) = -1.0  ! section averaged along z
987             ELSE
988                netcdf_data(i) = zu( section(i,1) )
989             ENDIF
990          ENDDO
991          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu_xy(av), &
992                                  netcdf_data, start = (/ 1 /),    &
993                                  count = (/ ns /) )
994          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 123 )
995
996!
997!--       Write zw data
998          DO  i = 1, ns
999             IF( section(i,1) == -1 )  THEN
1000                netcdf_data(i) = -1.0  ! section averaged along z
1001             ELSE
1002                netcdf_data(i) = zw( section(i,1) )
1003             ENDIF
1004          ENDDO
1005          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zw_xy(av), &
1006                                  netcdf_data, start = (/ 1 /),    &
1007                                  count = (/ ns /) )
1008          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 124 )
1009
1010!
1011!--       Write gridpoint number data
1012          netcdf_data(1:ns) = section(1:ns,1)
1013          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_ind_z_xy(av), &
1014                                  netcdf_data, start = (/ 1 /),       &
1015                                  count = (/ ns /) )
1016          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 125 )
1017
1018          DEALLOCATE( netcdf_data )
1019
1020!
1021!--       Write the cross section height u*, t*
1022          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu1_xy(av), &
1023                                  (/ zu(nzb+1) /), start = (/ 1 /), &
1024                                  count = (/ 1 /) )
1025          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 126 )
1026
1027!
1028!--       Write data for x and xu axis (shifted by -dx/2)
1029          ALLOCATE( netcdf_data(0:nx+1) )
1030
1031          DO  i = 0, nx+1
1032             netcdf_data(i) = i * dx
1033          ENDDO
1034
1035          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_x_xy(av), netcdf_data, &
1036                                  start = (/ 1 /), count = (/ nx+2 /) )
1037          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 127 )
1038
1039          DO  i = 0, nx+1
1040             netcdf_data(i) = ( i - 0.5 ) * dx
1041          ENDDO
1042
1043          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_xu_xy(av), &
1044                                  netcdf_data, start = (/ 1 /),    &
1045                                  count = (/ nx+2 /) )
1046          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 367 )
1047
1048          DEALLOCATE( netcdf_data )
1049
1050!
1051!--       Write data for y and yv axis (shifted by -dy/2)
1052          ALLOCATE( netcdf_data(0:ny+1) )
1053
1054          DO  i = 0, ny+1
1055             netcdf_data(i) = i * dy
1056          ENDDO
1057
1058          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_y_xy(av), netcdf_data, &
1059                                  start = (/ 1 /), count = (/ ny+2 /))
1060          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 128 )
1061
1062          DO  i = 0, ny+1
1063             netcdf_data(i) = ( i - 0.5 ) * dy
1064          ENDDO
1065
1066          nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_yv_xy(av), &
1067                                  netcdf_data, start = (/ 1 /),    &
1068                                  count = (/ ny+2 /))
1069          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 368 )
1070
1071          DEALLOCATE( netcdf_data )
1072
1073!
1074!--       In case of non-flat topography write height information
1075          IF ( TRIM( topography ) /= 'flat' )  THEN
1076
1077             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zusi_xy(av), &
1078                                     zu_s_inner(0:nx+1,0:ny+1), &
1079                                     start = (/ 1, 1 /), &
1080                                     count = (/ nx+2, ny+2 /) )
1081             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 427 )
1082
1083             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av), &
1084                                     zw_w_inner(0:nx+1,0:ny+1), &
1085                                     start = (/ 1, 1 /), &
1086                                     count = (/ nx+2, ny+2 /) )
1087             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 428 )
1088
1089          ENDIF
1090
1091
1092       CASE ( 'xy_ext' )
1093
1094!
1095!--       Get the list of variables and compare with the actual run.
1096!--       First var_list_old has to be reset, since GET_ATT does not assign
1097!--       trailing blanks.
1098          var_list_old = ' '
1099          nc_stat = NF90_GET_ATT( id_set_xy(av), NF90_GLOBAL, 'VAR_LIST', &
1100                                  var_list_old )
1101          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 129 )
1102
1103          var_list = ';'
1104          i = 1
1105          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
1106             IF ( INDEX( do2d(av,i), 'xy' ) /= 0 )  THEN
1107                netcdf_var_name = do2d(av,i)
1108                CALL clean_netcdf_varname( netcdf_var_name )
1109                var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
1110             ENDIF
1111             i = i + 1
1112          ENDDO
1113
1114          IF ( av == 0 )  THEN
1115             var = '(xy)'
1116          ELSE
1117             var = '(xy_av)'
1118          ENDIF
1119
1120          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
1121             PRINT*, '+++ WARNING: NetCDF file for cross-sections ' // &
1122                                   TRIM( var ) // ' from previuos run found,'
1123             PRINT*, '             but this file cannot be extended due to' // &
1124                                   ' variable mismatch.'
1125             PRINT*, '             New file is created instead.'
1126             extend = .FALSE.
1127             RETURN
1128          ENDIF
1129
1130!
1131!--       Calculate the number of current sections
1132          ns = 1
1133          DO WHILE ( section(ns,1) /= -9999  .AND.  ns <= 100 )
1134             ns = ns + 1
1135          ENDDO
1136          ns = ns - 1
1137
1138!
1139!--       Get and compare the number of horizontal cross sections
1140          nc_stat = NF90_INQ_VARID( id_set_xy(av), 'zu_xy', id_var_zu_xy(av) )
1141          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 130 )
1142
1143          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_zu_xy(av), &
1144                                           dimids = id_dim_zu_xy_old )
1145          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 131 )
1146          id_dim_zu_xy(av) = id_dim_zu_xy_old(1)
1147
1148          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_zu_xy(av), &
1149                                            len = ns_old )
1150          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 132 )
1151
1152          IF ( ns /= ns_old )  THEN
1153             PRINT*, '+++ WARNING: NetCDF file for cross-sections ' // &
1154                                   TRIM( var ) // ' from previuos run found,'
1155             PRINT*, '             but this file cannot be extended due to' // &
1156                                   ' mismatch in number of'
1157             PRINT*, '             cross sections.'
1158             PRINT*, '             New file is created instead.'
1159             extend = .FALSE.
1160             RETURN
1161          ENDIF
1162
1163!
1164!--       Get and compare the heights of the cross sections
1165          ALLOCATE( netcdf_data(1:ns_old) )
1166
1167          nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_zu_xy(av), netcdf_data )
1168          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 133 )
1169
1170          DO  i = 1, ns
1171             IF ( section(i,1) /= -1 )  THEN
1172                IF ( zu(section(i,1)) /= netcdf_data(i) )  THEN
1173                   PRINT*, '+++ WARNING: NetCDF file for cross-sections ' // &
1174                                     TRIM( var ) // ' from previuos run found,'
1175                   PRINT*, '             but this file cannot be extended' // &
1176                                      ' due to mismatch in cross'
1177                   PRINT*, '             section levels.'
1178                   PRINT*, '             New file is created instead.'
1179                   extend = .FALSE.
1180                   RETURN
1181                ENDIF
1182             ELSE
1183                IF ( -1.0 /= netcdf_data(i) )  THEN
1184                   PRINT*, '+++ WARNING: NetCDF file for cross-sections ' // &
1185                                     TRIM( var ) // ' from previuos run found,'
1186                   PRINT*, '             but this file cannot be extended' // &
1187                                      ' due to mismatch in cross'
1188                   PRINT*, '             section levels.'
1189                   PRINT*, '             New file is created instead.'
1190                   extend = .FALSE.
1191                   RETURN
1192                ENDIF
1193             ENDIF
1194          ENDDO
1195
1196          DEALLOCATE( netcdf_data )
1197
1198!
1199!--       Get the id of the time coordinate (unlimited coordinate) and its
1200!--       last index on the file. The next time level is do2d..count+1.
1201!--       The current time must be larger than the last output time
1202!--       on the file.
1203          nc_stat = NF90_INQ_VARID( id_set_xy(av), 'time', id_var_time_xy(av) )
1204          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 134 )
1205
1206          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_time_xy(av), &
1207                                           dimids = id_dim_time_old )
1208          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 135 )
1209          id_dim_time_xy(av) = id_dim_time_old(1)
1210
1211          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_time_xy(av), &
1212                                            len = do2d_xy_time_count(av) )
1213          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 136 )
1214
1215          nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_time_xy(av),    &
1216                                  last_time_coordinate,                 &
1217                                  start = (/ do2d_xy_time_count(av) /), &
1218                                  count = (/ 1 /) )
1219          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 137 )
1220
1221          IF ( last_time_coordinate(1) >= simulated_time )  THEN
1222             PRINT*, '+++ WARNING: NetCDF file for cross sections ' // &
1223                                   TRIM( var ) // ' from previuos run found,'
1224             PRINT*, '             but this file cannot be extended becaus' // &
1225                                   'e the current output time'
1226             PRINT*, '             is less or equal than the last output t' // &
1227                                   'ime on this file.'
1228             PRINT*, '             New file is created instead.'
1229             do2d_xy_time_count(av) = 0
1230             extend = .FALSE.
1231             RETURN
1232          ENDIF
1233
1234!
1235!--       Dataset seems to be extendable.
1236!--       Now get the variable ids.
1237          i = 1
1238          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
1239             IF ( INDEX( do2d(av,i), 'xy' ) /= 0 )  THEN
1240                netcdf_var_name = do2d(av,i)
1241                CALL clean_netcdf_varname( netcdf_var_name )
1242                nc_stat = NF90_INQ_VARID( id_set_xy(av), netcdf_var_name, &
1243                                          id_var_do2d(av,i) )
1244                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 138 )
1245             ENDIF
1246             i = i + 1
1247          ENDDO
1248
1249!
1250!--       Change the titel attribute on file
1251          nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'title', &
1252                                  TRIM( run_description_header ) )
1253          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 139 )
1254          PRINT*, '*** NetCDF file for cross-sections ' // TRIM( var ) // &
1255                       ' from previous run found.'
1256          PRINT*, '    This file will be extended.'
1257
1258
1259       CASE ( 'xz_new' )
1260
1261!
1262!--       Define some global attributes of the dataset
1263          nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'Conventions', &
1264                                  'COARDS' )
1265          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 140 )
1266
1267          IF ( av == 0 )  THEN
1268             time_average_text = ' '
1269          ELSE
1270             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
1271                                                            averaging_interval
1272          ENDIF
1273          nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'title', &
1274                                  TRIM( run_description_header )  //   &
1275                                  TRIM( time_average_text ) )
1276          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 141 )
1277          IF ( av == 1 )  THEN
1278             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
1279             nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'time_avg', &
1280                                     TRIM( time_average_text ) )
1281             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 141 )
1282          ENDIF
1283
1284!
1285!--       Define time coordinate for xz sections (unlimited dimension)
1286          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'time', NF90_UNLIMITED, &
1287                                  id_dim_time_xz(av) )
1288          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 142 )
1289
1290          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'time', NF90_DOUBLE, &
1291                                  id_dim_time_xz(av), id_var_time_xz(av) )
1292          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 143 )
1293
1294          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_time_xz(av), 'units', &
1295                                  'seconds')
1296          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 144 )
1297
1298!
1299!--       Define the spatial dimensions and coordinates for xz-sections.
1300!--       First, determine the number of vertical sections to be written.
1301          IF ( section(1,2) == -9999 )  THEN
1302             RETURN
1303          ELSE
1304             ns = 1
1305             DO WHILE ( section(ns,2) /= -9999  .AND.  ns <= 100 )
1306                ns = ns + 1
1307             ENDDO
1308             ns = ns - 1
1309          ENDIF
1310
1311!
1312!--       Define y-axis (for scalar position)
1313          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'y_xz', ns, id_dim_y_xz(av) )
1314          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 145 )
1315
1316          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'y_xz', NF90_DOUBLE, &
1317                                  id_dim_y_xz(av), id_var_y_xz(av) )
1318          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 146 )
1319
1320          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_y_xz(av), 'units', &
1321                                  'meters' )
1322          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 147 )
1323
1324!
1325!--       Define y-axis (for v position)
1326          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'yv_xz', ns, id_dim_yv_xz(av) )
1327          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 369 )
1328
1329          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'yv_xz', NF90_DOUBLE, &
1330                                  id_dim_yv_xz(av), id_var_yv_xz(av) )
1331          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 370 )
1332
1333          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_yv_xz(av), 'units', &
1334                                  'meters' )
1335          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 371 )
1336
1337!
1338!--       Define a variable to store the layer indices of the vertical cross
1339!--       sections
1340          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'ind_y_xz', NF90_DOUBLE, &
1341                                  id_dim_y_xz(av), id_var_ind_y_xz(av) )
1342          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 148 )
1343
1344          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_ind_y_xz(av), 'units', &
1345                                  'gridpoints')
1346          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 149 )
1347
1348!
1349!--       Define x-axis (for scalar position)
1350          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'x', nx+2, id_dim_x_xz(av) )
1351          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 150 )
1352
1353          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'x', NF90_DOUBLE, &
1354                                  id_dim_x_xz(av), id_var_x_xz(av) )
1355          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 151 )
1356
1357          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_x_xz(av), 'units', &
1358                                  'meters' )
1359          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 152 )
1360
1361!
1362!--       Define x-axis (for u position)
1363          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'xu', nx+2, id_dim_xu_xz(av) )
1364          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 372 )
1365
1366          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'xu', NF90_DOUBLE, &
1367                                  id_dim_xu_xz(av), id_var_xu_xz(av) )
1368          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 373 )
1369
1370          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_xu_xz(av), 'units', &
1371                                  'meters' )
1372          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 374 )
1373
1374!
1375!--       Define the two z-axes (zu and zw)
1376          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'zu', nz+2, id_dim_zu_xz(av) )
1377          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 153 )
1378
1379          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'zu', NF90_DOUBLE, &
1380                                  id_dim_zu_xz(av), id_var_zu_xz(av) )
1381          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 154 )
1382
1383          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_zu_xz(av), 'units', &
1384                                  'meters' )
1385          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 155 )
1386
1387          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'zw', nz+2, id_dim_zw_xz(av) )
1388          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 156 )
1389
1390          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'zw', NF90_DOUBLE, &
1391                                  id_dim_zw_xz(av), id_var_zw_xz(av) )
1392          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 157 )
1393
1394          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_zw_xz(av), 'units', &
1395                                  'meters' )
1396          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 158 )
1397
1398
1399!
1400!--       Define the variables
1401          var_list = ';'
1402          i = 1
1403
1404          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
1405
1406             IF ( INDEX( do2d(av,i), 'xz' ) /= 0 )  THEN
1407
1408!
1409!--             Check for the grid
1410                found = .TRUE.
1411                SELECT CASE ( do2d(av,i) )
1412!
1413!--                Most variables are defined on the zu grid
1414                   CASE ( 'e_xz', 'p_xz', 'pc_xz', 'pr_xz', 'pt_xz', 'q_xz',  &
1415                          'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qv_xz', &
1416                          's_xz', 'vpt_xz' )
1417
1418                      grid_x = 'x'
1419                      grid_y = 'y'
1420                      grid_z = 'zu'
1421!
1422!--                u grid
1423                   CASE ( 'u_xz' )
1424
1425                      grid_x = 'xu'
1426                      grid_y = 'y'
1427                      grid_z = 'zu'
1428!
1429!--                v grid
1430                   CASE ( 'v_xz' )
1431
1432                      grid_x = 'x'
1433                      grid_y = 'yv'
1434                      grid_z = 'zu'
1435!
1436!--                w grid
1437                   CASE ( 'w_xz' )
1438
1439                      grid_x = 'x'
1440                      grid_y = 'y'
1441                      grid_z = 'zw'
1442
1443                   CASE DEFAULT
1444!
1445!--                   Check for user-defined quantities
1446                      CALL user_define_netcdf_grid( do2d(av,i), found, &
1447                                                    grid_x, grid_y, grid_z )
1448
1449                END SELECT
1450
1451!
1452!--             Select the respective dimension ids
1453                IF ( grid_x == 'x' )  THEN
1454                   id_x = id_dim_x_xz(av)
1455                ELSEIF ( grid_x == 'xu' )  THEN
1456                   id_x = id_dim_xu_xz(av)
1457                ENDIF
1458
1459                IF ( grid_y == 'y' )  THEN
1460                   id_y = id_dim_y_xz(av)
1461                ELSEIF ( grid_y == 'yv' )  THEN
1462                   id_y = id_dim_yv_xz(av)
1463                ENDIF
1464
1465                IF ( grid_z == 'zu' )  THEN
1466                   id_z = id_dim_zu_xz(av)
1467                ELSEIF ( grid_z == 'zw' )  THEN
1468                   id_z = id_dim_zw_xz(av)
1469                ENDIF
1470
1471!
1472!--             Define the grid
1473                nc_stat = NF90_DEF_VAR( id_set_xz(av), do2d(av,i),             &
1474                                        nc_precision(2),                       &
1475                                   (/ id_x, id_y, id_z, id_dim_time_xz(av) /), &
1476                                        id_var_do2d(av,i) )
1477
1478                IF ( .NOT. found )  THEN
1479                   PRINT*, '+++ define_netcdf_header: no grid defined for', &
1480                                ' variable ', do2d(av,i)
1481                ENDIF
1482
1483                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
1484
1485                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 159 )
1486!
1487!--             Store the 'real' name of the variable (with *, for example)
1488!--             in the long_name attribute. This is evaluated by Ferret,
1489!--             for example.
1490                nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_do2d(av,i), &
1491                                        'long_name', do2d(av,i) )
1492                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 160 )
1493!
1494!--             Define the variable's unit
1495                nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_do2d(av,i), &
1496                                        'units', TRIM( do2d_unit(av,i) ) )
1497                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 355 )
1498             ENDIF
1499
1500             i = i + 1
1501
1502          ENDDO
1503
1504!
1505!--       No arrays to output. Close the netcdf file and return.
1506          IF ( i == 1 )  RETURN
1507
1508!
1509!--       Write the list of variables as global attribute (this is used by
1510!--       restart runs and by combine_plot_fields)
1511          nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'VAR_LIST', &
1512                                  var_list )
1513          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 161 )
1514
1515!
1516!--       Leave NetCDF define mode
1517          nc_stat = NF90_ENDDEF( id_set_xz(av) )
1518          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 162 )
1519
1520!
1521!--       Write axis data: y_xz, x, zu, zw
1522          ALLOCATE( netcdf_data(1:ns) )
1523
1524!
1525!--       Write y_xz data
1526          DO  i = 1, ns
1527             IF( section(i,2) == -1 )  THEN
1528                netcdf_data(i) = -1.0  ! section averaged along y
1529             ELSE
1530                netcdf_data(i) = section(i,2) * dy
1531             ENDIF
1532          ENDDO
1533          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_y_xz(av), netcdf_data, &
1534                                  start = (/ 1 /), count = (/ ns /) )
1535          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 163 )
1536
1537!
1538!--       Write yv_xz data
1539          DO  i = 1, ns
1540             IF( section(i,2) == -1 )  THEN
1541                netcdf_data(i) = -1.0  ! section averaged along y
1542             ELSE
1543                netcdf_data(i) = ( section(i,2) - 0.5 ) * dy
1544             ENDIF
1545          ENDDO
1546          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_yv_xz(av), &
1547                                  netcdf_data, start = (/ 1 /),    &
1548                                  count = (/ ns /) )
1549          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 375 )
1550
1551!
1552!--       Write gridpoint number data
1553          netcdf_data(1:ns) = section(1:ns,2)
1554          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_ind_y_xz(av), &
1555                                  netcdf_data, start = (/ 1 /),       &
1556                                  count = (/ ns /) )
1557          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 164 )
1558
1559          DEALLOCATE( netcdf_data )
1560
1561!
1562!--       Write data for x and xu axis (shifted by -dx/2)
1563          ALLOCATE( netcdf_data(0:nx+1) )
1564
1565          DO  i = 0, nx+1
1566             netcdf_data(i) = i * dx
1567          ENDDO
1568
1569          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_x_xz(av), netcdf_data, &
1570                                  start = (/ 1 /), count = (/ nx+2 /) )
1571          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 165 )
1572
1573          DO  i = 0, nx+1
1574             netcdf_data(i) = ( i - 0.5 ) * dx
1575          ENDDO
1576
1577          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_xu_xz(av), &
1578                                  netcdf_data, start = (/ 1 /),    &
1579                                  count = (/ nx+2 /) )
1580          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 376 )
1581
1582          DEALLOCATE( netcdf_data )
1583
1584!
1585!--       Write zu and zw data (vertical axes)
1586          ALLOCATE( netcdf_data(0:nz+1) )
1587
1588          netcdf_data(0:nz+1) = zu(nzb:nzt+1)
1589          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zu_xz(av), &
1590                                  netcdf_data, start = (/ 1 /),    &
1591                                  count = (/ nz+2 /) )
1592          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 166 )
1593
1594          netcdf_data(0:nz+1) = zw(nzb:nzt+1)
1595          nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zw_xz(av), &
1596                                  netcdf_data, start = (/ 1 /),    &
1597                                  count = (/ nz+2 /) )
1598          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 167 )
1599
1600          DEALLOCATE( netcdf_data )
1601
1602
1603       CASE ( 'xz_ext' )
1604
1605!
1606!--       Get the list of variables and compare with the actual run.
1607!--       First var_list_old has to be reset, since GET_ATT does not assign
1608!--       trailing blanks.
1609          var_list_old = ' '
1610          nc_stat = NF90_GET_ATT( id_set_xz(av), NF90_GLOBAL, 'VAR_LIST', &
1611                                  var_list_old )
1612          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 168 )
1613
1614          var_list = ';'
1615          i = 1
1616          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
1617             IF ( INDEX( do2d(av,i), 'xz' ) /= 0 )  THEN
1618                netcdf_var_name = do2d(av,i)
1619                CALL clean_netcdf_varname( netcdf_var_name )
1620                var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
1621             ENDIF
1622             i = i + 1
1623          ENDDO
1624
1625          IF ( av == 0 )  THEN
1626             var = '(xz)'
1627          ELSE
1628             var = '(xz_av)'
1629          ENDIF
1630
1631          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
1632             PRINT*, '+++ WARNING: NetCDF file for cross-sections ' // &
1633                                   TRIM( var ) // ' from previuos run found,'
1634             PRINT*, '             but this file cannot be extended due to' // &
1635                                   ' variable mismatch.'
1636             PRINT*, '             New file is created instead.'
1637             extend = .FALSE.
1638             RETURN
1639          ENDIF
1640
1641!
1642!--       Calculate the number of current sections
1643          ns = 1
1644          DO WHILE ( section(ns,2) /= -9999  .AND.  ns <= 100 )
1645             ns = ns + 1
1646          ENDDO
1647          ns = ns - 1
1648
1649!
1650!--       Get and compare the number of vertical cross sections
1651          nc_stat = NF90_INQ_VARID( id_set_xz(av), 'y_xz', id_var_y_xz(av) )
1652          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 169 )
1653
1654          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_y_xz(av), &
1655                                           dimids = id_dim_y_xz_old )
1656          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 170 )
1657          id_dim_y_xz(av) = id_dim_y_xz_old(1)
1658
1659          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_y_xz(av), &
1660                                            len = ns_old )
1661          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 171 )
1662
1663          IF ( ns /= ns_old )  THEN
1664             PRINT*, '+++ WARNING: NetCDF file for cross-sections ' // &
1665                                   TRIM( var ) // ' from previuos run found,'
1666             PRINT*, '             but this file cannot be extended due to' // &
1667                                   ' mismatch in number of'
1668             PRINT*, '             cross sections.'
1669             PRINT*, '             New file is created instead.'
1670             extend = .FALSE.
1671             RETURN
1672          ENDIF
1673
1674!
1675!--       Get and compare the heights of the cross sections
1676          ALLOCATE( netcdf_data(1:ns_old) )
1677
1678          nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_y_xz(av), netcdf_data )
1679          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 172 )
1680
1681          DO  i = 1, ns
1682             IF ( section(i,2) /= -1 )  THEN
1683                IF ( ( section(i,2) * dy ) /= netcdf_data(i) )  THEN
1684                   PRINT*, '+++ WARNING: NetCDF file for cross-sections ' // &
1685                                     TRIM( var ) // ' from previuos run found,'
1686                   PRINT*, '             but this file cannot be extended' // &
1687                                      ' due to mismatch in cross'
1688                   PRINT*, '             section indices.'
1689                   PRINT*, '             New file is created instead.'
1690                   extend = .FALSE.
1691                   RETURN
1692                ENDIF
1693             ELSE
1694                IF ( -1.0 /= netcdf_data(i) )  THEN
1695                   PRINT*, '+++ WARNING: NetCDF file for cross-sections ' // &
1696                                     TRIM( var ) // ' from previuos run found,'
1697                   PRINT*, '             but this file cannot be extended' // &
1698                                      ' due to mismatch in cross'
1699                   PRINT*, '             section indices.'
1700                   PRINT*, '             New file is created instead.'
1701                   extend = .FALSE.
1702                   RETURN
1703                ENDIF
1704             ENDIF
1705          ENDDO
1706
1707          DEALLOCATE( netcdf_data )
1708
1709!
1710!--       Get the id of the time coordinate (unlimited coordinate) and its
1711!--       last index on the file. The next time level is do2d..count+1.
1712!--       The current time must be larger than the last output time
1713!--       on the file.
1714          nc_stat = NF90_INQ_VARID( id_set_xz(av), 'time', id_var_time_xz(av) )
1715          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 173 )
1716
1717          nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_time_xz(av), &
1718                                           dimids = id_dim_time_old )
1719          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 174 )
1720          id_dim_time_xz(av) = id_dim_time_old(1)
1721
1722          nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_time_xz(av), &
1723                                            len = do2d_xz_time_count(av) )
1724          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 175 )
1725
1726          nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_time_xz(av),    &
1727                                  last_time_coordinate,                 &
1728                                  start = (/ do2d_xz_time_count(av) /), &
1729                                  count = (/ 1 /) )
1730          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 176 )
1731
1732          IF ( last_time_coordinate(1) >= simulated_time )  THEN
1733             PRINT*, '+++ WARNING: NetCDF file for cross sections ' // &
1734                                   TRIM( var ) // ' from previuos run found,'
1735             PRINT*, '             but this file cannot be extended becaus' // &
1736                                   'e the current output time'
1737             PRINT*, '             is less or equal than the last output t' // &
1738                                   'ime on this file.'
1739             PRINT*, '             New file is created instead.'
1740             do2d_xz_time_count(av) = 0
1741             extend = .FALSE.
1742             RETURN
1743          ENDIF
1744
1745!
1746!--       Dataset seems to be extendable.
1747!--       Now get the variable ids.
1748          i = 1
1749          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
1750             IF ( INDEX( do2d(av,i), 'xz' ) /= 0 )  THEN
1751                netcdf_var_name = do2d(av,i)
1752                CALL clean_netcdf_varname( netcdf_var_name )
1753                nc_stat = NF90_INQ_VARID( id_set_xz(av), netcdf_var_name, &
1754                                          id_var_do2d(av,i) )
1755                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 177 )
1756             ENDIF
1757             i = i + 1
1758          ENDDO
1759
1760!
1761!--       Change the titel attribute on file
1762          nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'title', &
1763                                  TRIM( run_description_header ) )
1764          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 178 )
1765          PRINT*, '*** NetCDF file for cross-sections ' // TRIM( var ) // &
1766                       ' from previous run found.'
1767          PRINT*, '    This file will be extended.'
1768
1769
1770       CASE ( 'yz_new' )
1771
1772!
1773!--       Define some global attributes of the dataset
1774          nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'Conventions', &
1775                                  'COARDS' )
1776          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 179 )
1777
1778          IF ( av == 0 )  THEN
1779             time_average_text = ' '
1780          ELSE
1781             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
1782                                                            averaging_interval
1783          ENDIF
1784          nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'title', &
1785                                  TRIM( run_description_header ) //    &
1786                                  TRIM( time_average_text ) )
1787          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 180 )
1788          IF ( av == 1 )  THEN
1789             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
1790             nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'time_avg', &
1791                                     TRIM( time_average_text ) )
1792             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 180 )
1793          ENDIF
1794
1795!
1796!--       Define time coordinate for yz sections (unlimited dimension)
1797          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'time', NF90_UNLIMITED, &
1798                                  id_dim_time_yz(av) )
1799          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 181 )
1800
1801          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'time', NF90_DOUBLE, &
1802                                  id_dim_time_yz(av), id_var_time_yz(av) )
1803          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 182 )
1804
1805          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_time_yz(av), 'units', &
1806                                  'seconds')
1807          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 183 )
1808
1809!
1810!--       Define the spatial dimensions and coordinates for yz-sections.
1811!--       First, determine the number of vertical sections to be written.
1812          IF ( section(1,3) == -9999 )  THEN
1813             RETURN
1814          ELSE
1815             ns = 1
1816             DO WHILE ( section(ns,3) /= -9999  .AND.  ns <= 100 )
1817                ns = ns + 1
1818             ENDDO
1819             ns = ns - 1
1820          ENDIF
1821
1822!
1823!--       Define x axis (for scalar position)
1824          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'x_yz', ns, id_dim_x_yz(av) )
1825          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 184 )
1826
1827          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'x_yz', NF90_DOUBLE, &
1828                                  id_dim_x_yz(av), id_var_x_yz(av) )
1829          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 185 )
1830
1831          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_x_yz(av), 'units', &
1832                                  'meters' )
1833          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 186 )
1834
1835!
1836!--       Define x axis (for u position)
1837          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'xu_yz', ns, id_dim_xu_yz(av) )
1838          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 377 )
1839
1840          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'xu_yz', NF90_DOUBLE, &
1841                                  id_dim_xu_yz(av), id_var_xu_yz(av) )
1842          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 378 )
1843
1844          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_xu_yz(av), 'units', &
1845                                  'meters' )
1846          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 379 )
1847
1848!
1849!--       Define a variable to store the layer indices of the vertical cross
1850!--       sections
1851          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'ind_x_yz', NF90_DOUBLE, &
1852                                  id_dim_x_yz(av), id_var_ind_x_yz(av) )
1853          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 187 )
1854
1855          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_ind_x_yz(av), 'units', &
1856                                  'gridpoints')
1857          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 188 )
1858
1859!
1860!--       Define y-axis (for scalar position)
1861          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'y', ny+2, id_dim_y_yz(av) )
1862          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 189 )
1863
1864          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'y', NF90_DOUBLE, &
1865                                  id_dim_y_yz(av), id_var_y_yz(av) )
1866          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 190 )
1867
1868          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_y_yz(av), 'units', &
1869                                  'meters' )
1870          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 191 )
1871
1872!
1873!--       Define y-axis (for v position)
1874          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'yv', ny+2, id_dim_yv_yz(av) )
1875          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 380 )
1876
1877          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'yv', NF90_DOUBLE, &
1878                                  id_dim_yv_yz(av), id_var_yv_yz(av) )
1879          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 381 )
1880
1881          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_yv_yz(av), 'units', &
1882                                  'meters' )
1883          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 382 )
1884
1885!
1886!--       Define the two z-axes (zu and zw)
1887          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'zu', nz+2, id_dim_zu_yz(av) )
1888          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 192 )
1889
1890          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'zu', NF90_DOUBLE, &
1891                                  id_dim_zu_yz(av), id_var_zu_yz(av) )
1892          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 193 )
1893
1894          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_zu_yz(av), 'units', &
1895                                  'meters' )
1896          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 194 )
1897
1898          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'zw', nz+2, id_dim_zw_yz(av) )
1899          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 195 )
1900
1901          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'zw', NF90_DOUBLE, &
1902                                  id_dim_zw_yz(av), id_var_zw_yz(av) )
1903          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 196 )
1904
1905          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_zw_yz(av), 'units', &
1906                                  'meters' )
1907          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 197 )
1908
1909
1910!
1911!--       Define the variables
1912          var_list = ';'
1913          i = 1
1914
1915          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
1916
1917             IF ( INDEX( do2d(av,i), 'yz' ) /= 0 )  THEN
1918
1919!
1920!--             Check for the grid
1921                found = .TRUE.
1922                SELECT CASE ( do2d(av,i) )
1923!
1924!--                Most variables are defined on the zu grid
1925                   CASE ( 'e_yz', 'p_yz', 'pc_yz', 'pr_yz', 'pt_yz', 'q_yz',  &
1926                          'ql_yz', 'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qv_yz', &
1927                          's_yz', 'vpt_yz' )
1928
1929                      grid_x = 'x'
1930                      grid_y = 'y'
1931                      grid_z = 'zu'
1932!
1933!--                u grid
1934                   CASE ( 'u_yz' )
1935
1936                      grid_x = 'xu'
1937                      grid_y = 'y'
1938                      grid_z = 'zu'
1939!
1940!--                v grid
1941                   CASE ( 'v_yz' )
1942
1943                      grid_x = 'x'
1944                      grid_y = 'yv'
1945                      grid_z = 'zu'
1946!
1947!--                w grid
1948                   CASE ( 'w_yz' )
1949
1950                      grid_x = 'x'
1951                      grid_y = 'y'
1952                      grid_z = 'zw'
1953
1954                   CASE DEFAULT
1955!
1956!--                   Check for user-defined quantities
1957                      CALL user_define_netcdf_grid( do2d(av,i), found, &
1958                                                    grid_x, grid_y, grid_z )
1959
1960                END SELECT
1961
1962!
1963!--             Select the respective dimension ids
1964                IF ( grid_x == 'x' )  THEN
1965                   id_x = id_dim_x_yz(av)
1966                ELSEIF ( grid_x == 'xu' )  THEN
1967                   id_x = id_dim_xu_yz(av)
1968                ENDIF
1969
1970                IF ( grid_y == 'y' )  THEN
1971                   id_y = id_dim_y_yz(av)
1972                ELSEIF ( grid_y == 'yv' )  THEN
1973                   id_y = id_dim_yv_yz(av)
1974                ENDIF
1975
1976                IF ( grid_z == 'zu' )  THEN
1977                   id_z = id_dim_zu_yz(av)
1978                ELSEIF ( grid_z == 'zw' )  THEN
1979                   id_z = id_dim_zw_yz(av)
1980                ENDIF
1981
1982!
1983!--             Define the grid
1984                nc_stat = NF90_DEF_VAR( id_set_yz(av), do2d(av,i),             &
1985                                        nc_precision(3),                       &
1986                                   (/ id_x, id_y, id_z, id_dim_time_yz(av) /), &
1987                                        id_var_do2d(av,i) )
1988
1989                IF ( .NOT. found )  THEN
1990                   PRINT*, '+++ define_netcdf_header: no grid defined for', &
1991                                ' variable ', do2d(av,i)
1992                ENDIF
1993
1994                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
1995
1996                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 198 )
1997!
1998!--             Store the 'real' name of the variable (with *, for example)
1999!--             in the long_name attribute. This is evaluated by Ferret,
2000!--             for example.
2001                nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_do2d(av,i), &
2002                                        'long_name', do2d(av,i) )
2003                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 199 )
2004!
2005!--             Define the variable's unit
2006                nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_do2d(av,i), &
2007                                        'units', TRIM( do2d_unit(av,i) ) )
2008                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 356 )
2009             ENDIF
2010
2011             i = i + 1
2012
2013          ENDDO
2014
2015!
2016!--       No arrays to output. Close the netcdf file and return.
2017          IF ( i == 1 )  RETURN
2018
2019!
2020!--       Write the list of variables as global attribute (this is used by
2021!--       restart runs and by combine_plot_fields)
2022          nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'VAR_LIST', &
2023                                  var_list )
2024          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 200 )
2025
2026!
2027!--       Leave NetCDF define mode
2028          nc_stat = NF90_ENDDEF( id_set_yz(av) )
2029          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 201 )
2030
2031!
2032!--       Write axis data: x_yz, y, zu, zw
2033          ALLOCATE( netcdf_data(1:ns) )
2034
2035!
2036!--       Write x_yz data
2037          DO  i = 1, ns
2038             IF( section(i,3) == -1 )  THEN
2039                netcdf_data(i) = -1.0  ! section averaged along x
2040             ELSE
2041                netcdf_data(i) = section(i,3) * dx
2042             ENDIF
2043          ENDDO
2044          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_x_yz(av), netcdf_data, &
2045                                  start = (/ 1 /), count = (/ ns /) )
2046          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 202 )
2047
2048!
2049!--       Write x_yz data (xu grid)
2050          DO  i = 1, ns
2051             IF( section(i,3) == -1 )  THEN
2052                netcdf_data(i) = -1.0  ! section averaged along x
2053             ELSE
2054                netcdf_data(i) = (section(i,3)-0.5) * dx
2055             ENDIF
2056          ENDDO
2057          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_xu_yz(av), netcdf_data, &
2058                                  start = (/ 1 /), count = (/ ns /) )
2059          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 383 )
2060
2061!
2062!--       Write gridpoint number data
2063          netcdf_data(1:ns) = section(1:ns,3)
2064          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_ind_x_yz(av), &
2065                                  netcdf_data, start = (/ 1 /),       &
2066                                  count = (/ ns /) )
2067          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 203 )
2068
2069          DEALLOCATE( netcdf_data )
2070
2071!
2072!--       Write data for y and yv axis (shifted by -dy/2)
2073          ALLOCATE( netcdf_data(0:ny+1) )
2074
2075          DO  j = 0, ny+1
2076             netcdf_data(j) = j * dy
2077          ENDDO
2078
2079          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_y_yz(av), netcdf_data, &
2080                                  start = (/ 1 /), count = (/ ny+2 /) )
2081          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 204 )
2082
2083          DO  j = 0, ny+1
2084             netcdf_data(j) = ( j - 0.5 ) * dy
2085          ENDDO
2086
2087          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_yv_yz(av), &
2088                                  netcdf_data, start = (/ 1 /),    &
2089                                  count = (/ ny+2 /) )
2090          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 384 )
2091
2092          DEALLOCATE( netcdf_data )
2093
2094!
2095!--       Write zu and zw data (vertical axes)
2096          ALLOCATE( netcdf_data(0:nz+1) )
2097
2098          netcdf_data(0:nz+1) = zu(nzb:nzt+1)
2099          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zu_yz(av), &
2100                                  netcdf_data, start = (/ 1 /),    &
2101                                  count = (/ nz+2 /) )
2102          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 205 )
2103
2104          netcdf_data(0:nz+1) = zw(nzb:nzt+1)
2105          nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_zw_yz(av), &
2106                                  netcdf_data, start = (/ 1 /),    &
2107                                  count = (/ nz+2 /) )
2108          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 206 )
2109
2110          DEALLOCATE( netcdf_data )
2111
2112
2113       CASE ( 'yz_ext' )
2114
2115!
2116!--       Get the list of variables and compare with the actual run.
2117!--       First var_list_old has to be reset, since GET_ATT does not assign
2118!--       trailing blanks.
2119          var_list_old = ' '
2120          nc_stat = NF90_GET_ATT( id_set_yz(av), NF90_GLOBAL, 'VAR_LIST', &
2121                                  var_list_old )
2122          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 207 )
2123
2124          var_list = ';'
2125          i = 1
2126          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2127             IF ( INDEX( do2d(av,i), 'yz' ) /= 0 )  THEN
2128                netcdf_var_name = do2d(av,i)
2129                CALL clean_netcdf_varname( netcdf_var_name )
2130                var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
2131             ENDIF
2132             i = i + 1
2133          ENDDO
2134
2135          IF ( av == 0 )  THEN
2136             var = '(yz)'
2137          ELSE
2138             var = '(yz_av)'
2139          ENDIF
2140
2141          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
2142             PRINT*, '+++ WARNING: NetCDF file for cross-sections ' // &
2143                                   TRIM( var ) // ' from previuos run found,'
2144             PRINT*, '             but this file cannot be extended due to' // &
2145                                   ' variable mismatch.'
2146             PRINT*, '             New file is created instead.'
2147             extend = .FALSE.
2148             RETURN
2149          ENDIF
2150
2151!
2152!--       Calculate the number of current sections
2153          ns = 1
2154          DO WHILE ( section(ns,3) /= -9999  .AND.  ns <= 100 )
2155             ns = ns + 1
2156          ENDDO
2157          ns = ns - 1
2158
2159!
2160!--       Get and compare the number of vertical cross sections
2161          nc_stat = NF90_INQ_VARID( id_set_yz(av), 'x_yz', id_var_x_yz(av) )
2162          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 208 )
2163
2164          nc_stat = NF90_INQUIRE_VARIABLE( id_set_yz(av), id_var_x_yz(av), &
2165                                           dimids = id_dim_x_yz_old )
2166          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 209 )
2167          id_dim_x_yz(av) = id_dim_x_yz_old(1)
2168
2169          nc_stat = NF90_INQUIRE_DIMENSION( id_set_yz(av), id_dim_x_yz(av), &
2170                                            len = ns_old )
2171          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 210 )
2172
2173          IF ( ns /= ns_old )  THEN
2174             PRINT*, '+++ WARNING: NetCDF file for cross-sections ' // &
2175                                   TRIM( var ) // ' from previuos run found,'
2176             PRINT*, '             but this file cannot be extended due to' // &
2177                                   ' mismatch in number of'
2178             PRINT*, '             cross sections.'
2179             PRINT*, '             New file is created instead.'
2180             extend = .FALSE.
2181             RETURN
2182          ENDIF
2183
2184!
2185!--       Get and compare the heights of the cross sections
2186          ALLOCATE( netcdf_data(1:ns_old) )
2187
2188          nc_stat = NF90_GET_VAR( id_set_yz(av), id_var_x_yz(av), netcdf_data )
2189          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 211 )
2190
2191          DO  i = 1, ns
2192             IF ( section(i,3) /= -1 )  THEN
2193                IF ( ( section(i,3) * dx ) /= netcdf_data(i) )  THEN
2194                   PRINT*, '+++ WARNING: NetCDF file for cross-sections ' // &
2195                                     TRIM( var ) // ' from previuos run found,'
2196                   PRINT*, '             but this file cannot be extended' // &
2197                                      ' due to mismatch in cross'
2198                   PRINT*, '             section indices.'
2199                   PRINT*, '             New file is created instead.'
2200                   extend = .FALSE.
2201                   RETURN
2202                ENDIF
2203             ELSE
2204                IF ( -1.0 /= netcdf_data(i) )  THEN
2205                   PRINT*, '+++ WARNING: NetCDF file for cross-sections ' // &
2206                                     TRIM( var ) // ' from previuos run found,'
2207                   PRINT*, '             but this file cannot be extended' // &
2208                                      ' due to mismatch in cross'
2209                   PRINT*, '             section indices.'
2210                   PRINT*, '             New file is created instead.'
2211                   extend = .FALSE.
2212                   RETURN
2213                ENDIF
2214             ENDIF
2215          ENDDO
2216
2217          DEALLOCATE( netcdf_data )
2218
2219!
2220!--       Get the id of the time coordinate (unlimited coordinate) and its
2221!--       last index on the file. The next time level is pl2d..count+1.
2222!--       The current time must be larger than the last output time
2223!--       on the file.
2224          nc_stat = NF90_INQ_VARID( id_set_yz(av), 'time', id_var_time_yz(av) )
2225          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 212 )
2226
2227          nc_stat = NF90_INQUIRE_VARIABLE( id_set_yz(av), id_var_time_yz(av), &
2228                                           dimids = id_dim_time_old )
2229          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 213 )
2230          id_dim_time_yz(av) = id_dim_time_old(1)
2231
2232          nc_stat = NF90_INQUIRE_DIMENSION( id_set_yz(av), id_dim_time_yz(av), &
2233                                            len = do2d_yz_time_count(av) )
2234          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 214 )
2235
2236          nc_stat = NF90_GET_VAR( id_set_yz(av), id_var_time_yz(av),    &
2237                                  last_time_coordinate,                 &
2238                                  start = (/ do2d_yz_time_count(av) /), &
2239                                  count = (/ 1 /) )
2240          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 215 )
2241
2242          IF ( last_time_coordinate(1) >= simulated_time )  THEN
2243             PRINT*, '+++ WARNING: NetCDF file for cross sections ' // &
2244                                   TRIM( var ) // ' from previuos run found,'
2245             PRINT*, '             but this file cannot be extended becaus' // &
2246                                   'e the current output time'
2247             PRINT*, '             is less or equal than the last output t' // &
2248                                   'ime on this file.'
2249             PRINT*, '             New file is created instead.'
2250             do2d_yz_time_count(av) = 0
2251             extend = .FALSE.
2252             RETURN
2253          ENDIF
2254
2255!
2256!--       Dataset seems to be extendable.
2257!--       Now get the variable ids.
2258          i = 1
2259          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2260             IF ( INDEX( do2d(av,i), 'yz' ) /= 0 )  THEN
2261                netcdf_var_name = do2d(av,i)
2262                CALL clean_netcdf_varname( netcdf_var_name )
2263                nc_stat = NF90_INQ_VARID( id_set_yz(av), netcdf_var_name, &
2264                                          id_var_do2d(av,i) )
2265                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 216 )
2266             ENDIF
2267             i = i + 1
2268          ENDDO
2269
2270!
2271!--       Change the titel attribute on file
2272          nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'title', &
2273                                  TRIM( run_description_header ) )
2274          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 217 )
2275          PRINT*, '*** NetCDF file for cross-sections ' // TRIM( var ) // &
2276                       ' from previous run found.'
2277          PRINT*, '    This file will be extended.'
2278
2279
2280       CASE ( 'pr_new' )
2281
2282!
2283!--       Define some global attributes of the dataset
2284          IF ( averaging_interval_pr /= 0.0 )  THEN
2285             WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
2286                                                            averaging_interval_pr
2287             nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'title',   &
2288                                     TRIM( run_description_header ) //  &
2289                                     TRIM( time_average_text ) )
2290             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 218 )
2291
2292             WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval_pr
2293             nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'time_avg', &
2294                                     TRIM( time_average_text ) )
2295          ELSE
2296             nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'title', &
2297                                     TRIM( run_description_header ) )
2298          ENDIF
2299          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 219 )
2300
2301!
2302!--       Define time coordinate for profiles (unlimited dimension)
2303          nc_stat = NF90_DEF_DIM( id_set_pr, 'time', NF90_UNLIMITED, &
2304                                  id_dim_time_pr )
2305          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 220 )
2306
2307          nc_stat = NF90_DEF_VAR( id_set_pr, 'time', NF90_DOUBLE, &
2308                                  id_dim_time_pr, id_var_time_pr )
2309          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 221 )
2310
2311          nc_stat = NF90_PUT_ATT( id_set_pr, id_var_time_pr, 'units', 'seconds')
2312          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 222 )
2313
2314!
2315!--       Define the variables
2316          var_list = ';'
2317          DO  i = 1, dopr_n
2318!
2319!--          First, remove those characters not allowed by NetCDF
2320             netcdf_var_name = data_output_pr(i)
2321             CALL clean_netcdf_varname( netcdf_var_name )
2322
2323             IF ( statistic_regions == 0 )  THEN
2324
2325!
2326!--             Define the z-axes (each variable gets its own z-axis)
2327                nc_stat = NF90_DEF_DIM( id_set_pr, 'z'//TRIM(netcdf_var_name), &
2328                                        nzt+2-nzb, id_dim_z_pr(i,0) )
2329                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 223 )
2330
2331                nc_stat = NF90_DEF_VAR( id_set_pr, 'z'//TRIM(netcdf_var_name), &
2332                                        NF90_DOUBLE, id_dim_z_pr(i,0),         &
2333                                        id_var_z_pr(i,0) )
2334                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 224 )
2335
2336                nc_stat = NF90_PUT_ATT( id_set_pr, id_var_z_pr(i,0), 'units', &
2337                                        'meters' )
2338                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 225 )
2339
2340!
2341!--             Define the variable
2342                nc_stat = NF90_DEF_VAR( id_set_pr, netcdf_var_name,           &
2343                                        nc_precision(5), (/ id_dim_z_pr(i,0), &
2344                                        id_dim_time_pr /), id_var_dopr(i,0) )
2345                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 226 )
2346
2347                nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,0), &
2348                                        'long_name', TRIM( data_output_pr(i) ) )
2349                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 227 )
2350
2351                nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,0), &
2352                                        'units', TRIM( dopr_unit(i) ) )
2353                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 228 )
2354
2355                var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
2356
2357             ELSE
2358!
2359!--             If statistic regions are defined, add suffix _SR+#SR to the
2360!--             names
2361                DO  j = 0, statistic_regions
2362                   WRITE ( suffix, '(''_'',I1)' )  j
2363
2364!
2365!--                Define the z-axes (each variable gets it own z-axis)
2366                   nc_stat = NF90_DEF_DIM( id_set_pr,                          &
2367                                           'z'//TRIM(netcdf_var_name)//suffix, &
2368                                           nzt+2-nzb, id_dim_z_pr(i,j) )
2369                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 229 )
2370
2371                   nc_stat = NF90_DEF_VAR( id_set_pr,                          &
2372                                           'z'//TRIM(netcdf_var_name)//suffix, &
2373                                           nc_precision(5), id_dim_z_pr(i,j),  &
2374                                           id_var_z_pr(i,j) )
2375                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 230 )
2376
2377                   nc_stat = NF90_PUT_ATT( id_set_pr, id_var_z_pr(i,j), &
2378                                           'units', 'meters' )
2379                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 231 )
2380
2381!
2382!--                Define the variable
2383                   nc_stat = NF90_DEF_VAR( id_set_pr,                         &
2384                                           TRIM( netcdf_var_name ) // suffix, &
2385                                           nc_precision(5),                   &
2386                                           (/ id_dim_z_pr(i,j),               &
2387                                           id_dim_time_pr /), id_var_dopr(i,j) )
2388                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 232 )
2389
2390                   nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,j),        &
2391                                           'long_name',                        &
2392                                           TRIM( data_output_pr(i) ) // ' SR ' &
2393                                           // suffix(2:2) )
2394                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 233 )
2395
2396                   nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,j), &
2397                                           'units', TRIM( dopr_unit(i) ) )
2398                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 234 )
2399
2400                   var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
2401                              suffix // ';'
2402
2403                ENDDO
2404
2405             ENDIF
2406
2407          ENDDO
2408
2409!
2410!--       Write the list of variables as global attribute (this is used by
2411!--       restart runs)
2412          nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'VAR_LIST', var_list )
2413          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 235 )
2414
2415!
2416!--       Define normalization variables (as time series)
2417          DO  i = 1, dopr_norm_num
2418
2419             nc_stat = NF90_DEF_VAR( id_set_pr, 'NORM_' // &
2420                                     TRIM( dopr_norm_names(i) ), &
2421                                     nc_precision(5), (/ id_dim_time_pr /), &
2422                                     id_var_norm_dopr(i) )
2423             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 236 )
2424
2425             nc_stat = NF90_PUT_ATT( id_set_pr, id_var_norm_dopr(i), &
2426                                     'long_name',                    &
2427                                     TRIM( dopr_norm_longnames(i) ) )
2428             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 237 )
2429
2430          ENDDO
2431
2432!
2433!--       Leave NetCDF define mode
2434          nc_stat = NF90_ENDDEF( id_set_pr )
2435          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 238 )
2436
2437!
2438!--       Write z-axes data
2439          DO  i = 1, dopr_n
2440             DO  j = 0, statistic_regions
2441
2442                nc_stat = NF90_PUT_VAR( id_set_pr, id_var_z_pr(i,j),      &
2443                                        hom(nzb:nzt+1,2,dopr_index(i),0), &
2444                                        start = (/ 1 /),                  &
2445                                        count = (/ nzt-nzb+2 /) )
2446                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 239 )
2447
2448             ENDDO
2449          ENDDO
2450
2451
2452       CASE ( 'pr_ext' )
2453
2454!
2455!--       Get the list of variables and compare with the actual run.
2456!--       First var_list_old has to be reset, since GET_ATT does not assign
2457!--       trailing blanks.
2458          var_list_old = ' '
2459          nc_stat = NF90_GET_ATT( id_set_pr, NF90_GLOBAL, 'VAR_LIST', &
2460                                  var_list_old )
2461          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 240 )
2462
2463          var_list = ';'
2464          DO  i = 1, dopr_n
2465
2466             netcdf_var_name = data_output_pr(i)
2467             CALL clean_netcdf_varname( netcdf_var_name )
2468
2469             IF ( statistic_regions == 0 )  THEN
2470                var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
2471             ELSE
2472                DO  j = 0, statistic_regions
2473                   WRITE ( suffix, '(''_'',I1)' )  j
2474                   var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
2475                              suffix // ';'
2476                ENDDO
2477             ENDIF
2478
2479          ENDDO
2480
2481          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
2482             PRINT*, '+++ WARNING: NetCDF file for vertical profiles from' // &
2483                                   ' previuos run found,'
2484             PRINT*, '             but this file cannot be extended due to' // &
2485                                   ' variable mismatch.'
2486             PRINT*, '             New file is created instead.'
2487             extend = .FALSE.
2488             RETURN
2489          ENDIF
2490
2491!
2492!--       Get the id of the time coordinate (unlimited coordinate) and its
2493!--       last index on the file. The next time level is dopr..count+1.
2494!--       The current time must be larger than the last output time
2495!--       on the file.
2496          nc_stat = NF90_INQ_VARID( id_set_pr, 'time', id_var_time_pr )
2497          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 241 )
2498
2499          nc_stat = NF90_INQUIRE_VARIABLE( id_set_pr, id_var_time_pr, &
2500                                           dimids = id_dim_time_old )
2501          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 242 )
2502          id_dim_time_pr = id_dim_time_old(1)
2503
2504          nc_stat = NF90_INQUIRE_DIMENSION( id_set_pr, id_dim_time_pr, &
2505                                            len = dopr_time_count )
2506          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 243 )
2507
2508          nc_stat = NF90_GET_VAR( id_set_pr, id_var_time_pr,        &
2509                                  last_time_coordinate,             &
2510                                  start = (/ dopr_time_count /), &
2511                                  count = (/ 1 /) )
2512          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 244 )
2513
2514          IF ( last_time_coordinate(1) >= simulated_time )  THEN
2515             PRINT*, '+++ WARNING: NetCDF file for vertical profiles from' // &
2516                                   ' previuos run found,'
2517             PRINT*, '             but this file cannot be extended becaus' // &
2518                                   'e the current output time'
2519             PRINT*, '             is less or equal than the last output t' // &
2520                                   'ime on this file.'
2521             PRINT*, '             New file is created instead.'
2522             dopr_time_count = 0
2523             extend = .FALSE.
2524             RETURN
2525          ENDIF
2526
2527!
2528!--       Dataset seems to be extendable.
2529!--       Now get the variable ids.
2530          i = 1
2531          DO  i = 1, dopr_n
2532 
2533             netcdf_var_name_base = data_output_pr(i)
2534             CALL clean_netcdf_varname( netcdf_var_name_base )
2535
2536             IF ( statistic_regions == 0 )  THEN
2537                nc_stat = NF90_INQ_VARID( id_set_pr, netcdf_var_name_base, &
2538                                          id_var_dopr(i,0) )
2539                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 245 )
2540             ELSE
2541                DO  j = 0, statistic_regions
2542                   WRITE ( suffix, '(''_'',I1)' )  j
2543                   netcdf_var_name = TRIM( netcdf_var_name_base ) // suffix
2544                   nc_stat = NF90_INQ_VARID( id_set_pr, netcdf_var_name, &
2545                                             id_var_dopr(i,j) )
2546                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 246 )
2547                ENDDO
2548             ENDIF
2549
2550          ENDDO
2551
2552!
2553!--       Get ids of the normalization variables
2554          DO  i = 1, dopr_norm_num
2555             nc_stat = NF90_INQ_VARID( id_set_pr,                             &
2556                                       'NORM_' // TRIM( dopr_norm_names(i) ), &
2557                                       id_var_norm_dopr(i) )
2558             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 247 )
2559          ENDDO
2560
2561!
2562!--       Change the title attribute on file
2563          nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'title', &
2564                                  TRIM( run_description_header ) )
2565          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 248 )
2566          PRINT*, '*** NetCDF file for vertical profiles from previous run' // &
2567                       ' found.'
2568          PRINT*, '    This file will be extended.'
2569
2570
2571       CASE ( 'ts_new' )
2572
2573!
2574!--       Define some global attributes of the dataset
2575          nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'title', &
2576                                  TRIM( run_description_header ) )
2577          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 249 )
2578
2579!
2580!--       Define time coordinate for time series (unlimited dimension)
2581          nc_stat = NF90_DEF_DIM( id_set_ts, 'time', NF90_UNLIMITED, &
2582                                  id_dim_time_ts )
2583          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 250 )
2584
2585          nc_stat = NF90_DEF_VAR( id_set_ts, 'time', NF90_DOUBLE, &
2586                                  id_dim_time_ts, id_var_time_ts )
2587          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 251 )
2588
2589          nc_stat = NF90_PUT_ATT( id_set_ts, id_var_time_ts, 'units', 'seconds')
2590          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 252 )
2591
2592!
2593!--       Define the variables
2594          var_list = ';'
2595          DO  i = 1, dots_num
2596!
2597!--          First, remove those characters not allowed by NetCDF
2598             netcdf_var_name = dots_label(i)
2599             CALL clean_netcdf_varname( netcdf_var_name )
2600
2601             IF ( statistic_regions == 0 )  THEN
2602
2603                nc_stat = NF90_DEF_VAR( id_set_ts, netcdf_var_name,            &
2604                                        nc_precision(6), (/ id_dim_time_ts /), &
2605                                        id_var_dots(i,0) )
2606                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 253 )
2607
2608                nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,0), &
2609                                        'long_name', TRIM( dots_label(i) ) )
2610                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 254 )
2611
2612                nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,0), &
2613                                        'units', TRIM( dots_unit(i) ) )
2614                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 255 )
2615
2616                var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
2617
2618             ELSE
2619!
2620!--             If statistic regions are defined, add suffix _SR+#SR to the
2621!--             names
2622                DO  j = 0, statistic_regions
2623                   WRITE ( suffix, '(''_'',I1)' )  j
2624
2625                   nc_stat = NF90_DEF_VAR( id_set_ts,                         &
2626                                           TRIM( netcdf_var_name ) // suffix, &
2627                                           nc_precision(6),                   &
2628                                           (/ id_dim_time_ts /),              &
2629                                           id_var_dots(i,j) )
2630                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 256 )
2631
2632                   nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,j),       &
2633                                           'long_name',                       &
2634                                           TRIM( dots_label(i) ) // ' SR ' // &
2635                                           suffix(2:2) )
2636                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 257 )
2637
2638                   var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
2639                              suffix // ';'
2640
2641                ENDDO
2642
2643             ENDIF
2644
2645          ENDDO
2646
2647!
2648!--       Write the list of variables as global attribute (this is used by
2649!--       restart runs)
2650          nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'VAR_LIST', var_list )
2651          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 258 )
2652
2653!
2654!--       Leave NetCDF define mode
2655          nc_stat = NF90_ENDDEF( id_set_ts )
2656          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 259 )
2657
2658
2659       CASE ( 'ts_ext' )
2660
2661!
2662!--       Get the list of variables and compare with the actual run.
2663!--       First var_list_old has to be reset, since GET_ATT does not assign
2664!--       trailing blanks.
2665          var_list_old = ' '
2666          nc_stat = NF90_GET_ATT( id_set_ts, NF90_GLOBAL, 'VAR_LIST', &
2667                                  var_list_old )
2668          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 260 )
2669
2670          var_list = ';'
2671          i = 1
2672          DO  i = 1, dots_num
2673
2674             netcdf_var_name = dots_label(i)
2675             CALL clean_netcdf_varname( netcdf_var_name )
2676
2677             IF ( statistic_regions == 0 )  THEN
2678                var_list = TRIM(var_list) // TRIM(netcdf_var_name) // ';'
2679             ELSE
2680                DO  j = 0, statistic_regions
2681                   WRITE ( suffix, '(''_'',I1)' )  j
2682                   var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
2683                              suffix // ';'
2684                ENDDO
2685             ENDIF
2686
2687          ENDDO
2688
2689          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
2690             PRINT*, '+++ WARNING: NetCDF file for time series from' //&
2691                                   ' previuos run found,'
2692             PRINT*, '             but this file cannot be extended due to' // &
2693                                   ' variable mismatch.'
2694             PRINT*, '             New file is created instead.'
2695             extend = .FALSE.
2696             RETURN
2697          ENDIF
2698
2699!
2700!--       Get the id of the time coordinate (unlimited coordinate) and its
2701!--       last index on the file. The next time level is dots..count+1.
2702!--       The current time must be larger than the last output time
2703!--       on the file.
2704          nc_stat = NF90_INQ_VARID( id_set_ts, 'time', id_var_time_ts )
2705          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 261 )
2706
2707          nc_stat = NF90_INQUIRE_VARIABLE( id_set_ts, id_var_time_ts, &
2708                                           dimids = id_dim_time_old )
2709          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 262 )
2710          id_dim_time_ts = id_dim_time_old(1)
2711
2712          nc_stat = NF90_INQUIRE_DIMENSION( id_set_ts, id_dim_time_ts, &
2713                                            len = dots_time_count )
2714          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 263 )
2715
2716          nc_stat = NF90_GET_VAR( id_set_ts, id_var_time_ts,        &
2717                                  last_time_coordinate,             &
2718                                  start = (/ dots_time_count /), &
2719                                  count = (/ 1 /) )
2720          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 264 )
2721
2722          IF ( last_time_coordinate(1) >= simulated_time )  THEN
2723             PRINT*, '+++ WARNING: NetCDF file for time series from' // &
2724                                   ' previuos run found,'
2725             PRINT*, '             but this file cannot be extended becaus' // &
2726                                   'e the current output time'
2727             PRINT*, '             is less or equal than the last output t' // &
2728                                   'ime on this file.'
2729             PRINT*, '             New file is created instead.'
2730             dots_time_count = 0
2731             extend = .FALSE.
2732             RETURN
2733          ENDIF
2734
2735!
2736!--       Dataset seems to be extendable.
2737!--       Now get the variable ids
2738          i = 1
2739          DO  i = 1, dots_num
2740 
2741             netcdf_var_name_base = dots_label(i)
2742             CALL clean_netcdf_varname( netcdf_var_name_base )
2743
2744             IF ( statistic_regions == 0 )  THEN
2745                nc_stat = NF90_INQ_VARID( id_set_ts, netcdf_var_name_base, &
2746                                          id_var_dots(i,0) )
2747                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 265 )
2748             ELSE
2749                DO  j = 0, statistic_regions
2750                   WRITE ( suffix, '(''_'',I1)' )  j
2751                   netcdf_var_name = TRIM( netcdf_var_name_base ) // suffix
2752                   nc_stat = NF90_INQ_VARID( id_set_ts, netcdf_var_name, &
2753                                             id_var_dots(i,j) )
2754                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 266 )
2755                ENDDO
2756             ENDIF
2757
2758          ENDDO
2759
2760!
2761!--       Change the title attribute on file
2762          nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'title', &
2763                                  TRIM( run_description_header ) )
2764          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 267 )
2765          PRINT*, '*** NetCDF file for time series from previous run found.'
2766          PRINT*, '    This file will be extended.'
2767
2768
2769       CASE ( 'sp_new' )
2770
2771!
2772!--       Define some global attributes of the dataset
2773          IF ( averaging_interval_sp /= 0.0 )  THEN
2774             WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
2775                                                            averaging_interval_sp
2776             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title',  &
2777                                     TRIM( run_description_header ) // &
2778                                     TRIM( time_average_text ) )
2779             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 268 )
2780
2781             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval_sp
2782             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'time_avg', &
2783                                     TRIM( time_average_text ) )
2784          ELSE
2785             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
2786                                     TRIM( run_description_header ) )
2787          ENDIF
2788          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 269 )
2789
2790!
2791!--       Define time coordinate for spectra (unlimited dimension)
2792          nc_stat = NF90_DEF_DIM( id_set_sp, 'time', NF90_UNLIMITED, &
2793                                  id_dim_time_sp )
2794          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 270 )
2795
2796          nc_stat = NF90_DEF_VAR( id_set_sp, 'time', NF90_DOUBLE, &
2797                                  id_dim_time_sp, id_var_time_sp )
2798          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 271 )
2799
2800          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_time_sp, 'units', 'seconds')
2801          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 272 )
2802
2803!
2804!--       Define the spatial dimensions and coordinates for spectra.
2805!--       First, determine the number of vertical levels for which spectra
2806!--       are to be output.
2807          ns = 1
2808          DO WHILE ( comp_spectra_level(ns) /= 999999  .AND.  ns <= 100 )
2809             ns = ns + 1
2810          ENDDO
2811          ns = ns - 1
2812
2813!
2814!--       Define vertical coordinate grid (zu grid)
2815          nc_stat = NF90_DEF_DIM( id_set_sp, 'zu_sp', ns, id_dim_zu_sp )
2816          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 273 )
2817
2818          nc_stat = NF90_DEF_VAR( id_set_sp, 'zu_sp', NF90_DOUBLE, &
2819                                  id_dim_zu_sp, id_var_zu_sp )
2820          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 274 )
2821
2822          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_zu_sp, 'units', 'meters' )
2823          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 275 )
2824
2825!
2826!--       Define vertical coordinate grid (zw grid)
2827          nc_stat = NF90_DEF_DIM( id_set_sp, 'zw_sp', ns, id_dim_zw_sp )
2828          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 276 )
2829
2830          nc_stat = NF90_DEF_VAR( id_set_sp, 'zw_sp', NF90_DOUBLE, &
2831                                  id_dim_zw_sp, id_var_zw_sp )
2832          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 277 )
2833
2834          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_zw_sp, 'units', 'meters' )
2835          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 278 )
2836
2837!
2838!--       Define x-axis
2839          nc_stat = NF90_DEF_DIM( id_set_sp, 'k_x', nx/2, id_dim_x_sp )
2840          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 279 )
2841
2842          nc_stat = NF90_DEF_VAR( id_set_sp, 'k_x', NF90_DOUBLE, id_dim_x_sp, &
2843                                  id_var_x_sp )
2844          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 280 )
2845
2846          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_x_sp, 'units', 'm-1' )
2847          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 281 )
2848
2849!
2850!--       Define y-axis
2851          nc_stat = NF90_DEF_DIM( id_set_sp, 'k_y', ny/2, id_dim_y_sp )
2852          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 282 )
2853
2854          nc_stat = NF90_DEF_VAR( id_set_sp, 'k_y', NF90_DOUBLE, id_dim_y_sp, &
2855                                  id_var_y_sp )
2856          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 283 )
2857
2858          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_y_sp, 'units', 'm-1' )
2859          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 284 )
2860
2861!
2862!--       Define the variables
2863          var_list = ';'
2864          i = 1
2865          DO WHILE ( data_output_sp(i) /= ' '  .AND.  i <= 10 )
2866
2867             IF ( INDEX( spectra_direction(i), 'x' ) /= 0 )  THEN
2868
2869!
2870!--             Define the variable
2871                netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
2872                IF ( data_output_sp(i) == 'w' )  THEN
2873                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
2874                                           nc_precision(7), (/ id_dim_x_sp, &
2875                                           id_dim_zw_sp, id_dim_time_sp /), &
2876                                           id_var_dospx(i) )
2877                ELSE
2878                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
2879                                           nc_precision(7), (/ id_dim_x_sp, &
2880                                           id_dim_zu_sp, id_dim_time_sp /), &
2881                                           id_var_dospx(i) )
2882                ENDIF
2883                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 285 )
2884
2885                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospx(i), &
2886                                        'long_name', netcdf_var_name )
2887                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 286 )
2888
2889                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospx(i), &
2890                                        'units', 'unknown' )
2891                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 287 )
2892
2893                var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
2894
2895             ENDIF
2896
2897             IF ( INDEX( spectra_direction(i), 'y' ) /= 0 )  THEN
2898
2899!
2900!--             Define the variable
2901                netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
2902                IF ( data_output_sp(i) == 'w' )  THEN
2903                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
2904                                           nc_precision(7), (/ id_dim_y_sp, &
2905                                           id_dim_zw_sp, id_dim_time_sp /), &
2906                                           id_var_dospy(i) )
2907                ELSE
2908                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
2909                                           nc_precision(7), (/ id_dim_y_sp, &
2910                                           id_dim_zu_sp, id_dim_time_sp /), &
2911                                           id_var_dospy(i) )
2912                ENDIF
2913                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 288 )
2914
2915                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospy(i), &
2916                                        'long_name', netcdf_var_name )
2917                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 289 )
2918
2919                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospy(i), &
2920                                        'units', 'unknown' )
2921                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 290 )
2922
2923                var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
2924
2925             ENDIF
2926
2927             i = i + 1
2928
2929          ENDDO
2930
2931!
2932!--       Write the list of variables as global attribute (this is used by
2933!--       restart runs)
2934          nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', var_list )
2935          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 291 )
2936
2937!
2938!--       Leave NetCDF define mode
2939          nc_stat = NF90_ENDDEF( id_set_sp )
2940          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 292 )
2941
2942!
2943!--       Write axis data: zu_sp, zw_sp, k_x, k_y
2944          ALLOCATE( netcdf_data(1:ns) )
2945
2946!
2947!--       Write zu data
2948          netcdf_data(1:ns) = zu( comp_spectra_level(1:ns) )
2949          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zu_sp, netcdf_data, &
2950                                  start = (/ 1 /), count = (/ ns /) )
2951          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 293 )
2952
2953!
2954!--       Write zw data
2955          netcdf_data(1:ns) = zw( comp_spectra_level(1:ns) )
2956          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zw_sp, netcdf_data, &
2957                                  start = (/ 1 /), count = (/ ns /) )
2958          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 294 )
2959
2960          DEALLOCATE( netcdf_data )
2961
2962!
2963!--       Write data for x and y axis (wavenumbers)
2964          ALLOCATE( netcdf_data(nx/2) )
2965          DO  i = 1, nx/2
2966             netcdf_data(i) = 2.0 * pi * i / ( dx * ( nx + 1 ) )
2967          ENDDO
2968
2969          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_x_sp, netcdf_data, &
2970                                  start = (/ 1 /), count = (/ nx/2 /) )
2971          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 295 )
2972
2973          DEALLOCATE( netcdf_data )
2974
2975          ALLOCATE( netcdf_data(ny/2) )
2976          DO  i = 1, ny/2
2977             netcdf_data(i) = 2.0 * pi * i / ( dy * ( ny + 1 ) )
2978          ENDDO
2979
2980          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_y_sp, netcdf_data, &
2981                                  start = (/ 1 /), count = (/ ny/2 /) )
2982          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 296 )
2983
2984          DEALLOCATE( netcdf_data )
2985
2986
2987       CASE ( 'sp_ext' )
2988
2989!
2990!--       Get the list of variables and compare with the actual run.
2991!--       First var_list_old has to be reset, since GET_ATT does not assign
2992!--       trailing blanks.
2993          var_list_old = ' '
2994          nc_stat = NF90_GET_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', &
2995                                  var_list_old )
2996          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 297 )
2997
2998          var_list = ';'
2999          i = 1
3000          DO WHILE ( data_output_sp(i) /= ' '  .AND.  i <= 10 )
3001
3002             IF ( INDEX( spectra_direction(i), 'x' ) /= 0 )  THEN
3003                netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
3004                var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
3005             ENDIF
3006
3007             IF ( INDEX( spectra_direction(i), 'y' ) /= 0 )  THEN
3008                netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
3009                var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
3010             ENDIF
3011
3012             i = i + 1
3013
3014          ENDDO
3015
3016          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
3017             PRINT*, '+++ WARNING: NetCDF file for spectra from previous ' //&
3018                                   'run found,'
3019             PRINT*, '             but this file cannot be extended due to' // &
3020                                   ' variable mismatch.'
3021             PRINT*, '             New file is created instead.'
3022             extend = .FALSE.
3023             RETURN
3024          ENDIF
3025
3026!
3027!--       Determine the number of current vertical levels for which spectra
3028!--       shall be output
3029          ns = 1
3030          DO WHILE ( comp_spectra_level(ns) /= 999999  .AND.  ns <= 100 )
3031             ns = ns + 1
3032          ENDDO
3033          ns = ns - 1
3034
3035!
3036!--       Get and compare the number of vertical levels
3037          nc_stat = NF90_INQ_VARID( id_set_sp, 'zu_sp', id_var_zu_sp )
3038          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 298 )
3039
3040          nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_zu_sp, &
3041                                           dimids = id_dim_zu_sp_old )
3042          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 299 )
3043          id_dim_zu_sp = id_dim_zu_sp_old(1)
3044
3045          nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_zu_sp, &
3046                                            len = ns_old )
3047          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 300 )
3048
3049          IF ( ns /= ns_old )  THEN
3050             PRINT*, '+++ WARNING: NetCDF file for spectra from previous ' //&
3051                                   'run found,'
3052             PRINT*, '             but this file cannot be extended due to' // &
3053                                   ' mismatch in number of'
3054             PRINT*, '             vertical levels.'
3055             PRINT*, '             New file is created instead.'
3056             extend = .FALSE.
3057             RETURN
3058          ENDIF
3059
3060!
3061!--       Get and compare the heights of the cross sections
3062          ALLOCATE( netcdf_data(1:ns_old) )
3063
3064          nc_stat = NF90_GET_VAR( id_set_sp, id_var_zu_sp, netcdf_data )
3065          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 301 )
3066
3067          DO  i = 1, ns
3068             IF ( zu(comp_spectra_level(i)) /= netcdf_data(i) )  THEN
3069                PRINT*, '+++ WARNING: NetCDF file for spectra from previou' // &
3070                                   's run found,'
3071                PRINT*, '             but this file cannot be extended due' // &
3072                                   ' to mismatch in heights'
3073                PRINT*, '             of vertical levels.'
3074                PRINT*, '             New file is created instead.'
3075                extend = .FALSE.
3076                RETURN
3077             ENDIF
3078          ENDDO
3079
3080          DEALLOCATE( netcdf_data )
3081
3082!
3083!--       Get the id of the time coordinate (unlimited coordinate) and its
3084!--       last index on the file. The next time level is plsp..count+1.
3085!--       The current time must be larger than the last output time
3086!--       on the file.
3087          nc_stat = NF90_INQ_VARID( id_set_sp, 'time', id_var_time_sp )
3088          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 302 )
3089
3090          nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_time_sp, &
3091                                           dimids = id_dim_time_old )
3092          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 303 )
3093          id_dim_time_sp = id_dim_time_old(1)
3094
3095          nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_time_sp, &
3096                                            len = dosp_time_count )
3097          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 304 )
3098
3099          nc_stat = NF90_GET_VAR( id_set_sp, id_var_time_sp,        &
3100                                  last_time_coordinate,             &
3101                                  start = (/ dosp_time_count /), &
3102                                  count = (/ 1 /) )
3103          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 305 )
3104
3105          IF ( last_time_coordinate(1) >= simulated_time )  THEN
3106             PRINT*, '+++ WARNING: NetCDF file for spectra from previous ' // &
3107                                   'run found,'
3108             PRINT*, '             but this file cannot be extended becaus' // &
3109                                   'e the current output time'
3110             PRINT*, '             is less or equal than the last output t' // &
3111                                   'ime on this file.'
3112             PRINT*, '             New file is created instead.'
3113             dosp_time_count = 0
3114             extend = .FALSE.
3115             RETURN
3116          ENDIF
3117
3118!
3119!--       Dataset seems to be extendable.
3120!--       Now get the variable ids.
3121          i = 1
3122          DO WHILE ( data_output_sp(i) /= ' '  .AND.  i <= 10 )
3123
3124             IF ( INDEX( spectra_direction(i), 'x' ) /= 0 )  THEN
3125                netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
3126                nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
3127                                          id_var_dospx(i) )
3128                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 306 )
3129             ENDIF
3130
3131             IF ( INDEX( spectra_direction(i), 'y' ) /= 0 )  THEN
3132                netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
3133                nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
3134                                          id_var_dospy(i) )
3135                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 307 )
3136             ENDIF
3137
3138             i = i + 1
3139
3140          ENDDO
3141
3142!
3143!--       Change the titel attribute on file
3144          IF ( averaging_interval_sp /= 0.0 )  THEN
3145             WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
3146                                                            averaging_interval_sp
3147             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title',  &
3148                                     TRIM( run_description_header ) // &
3149                                     TRIM( time_average_text ) )
3150             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 308 )
3151
3152             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval_sp
3153             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'time_avg', &
3154                                     TRIM( time_average_text ) )
3155          ELSE
3156             nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'title', &
3157                                     TRIM( run_description_header ) )
3158          ENDIF
3159          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 309 )
3160
3161          PRINT*, '*** NetCDF file for spectra from previous run found.'
3162          PRINT*, '    This file will be extended.'
3163
3164
3165       CASE ( 'pt_new' )
3166
3167!
3168!--       Define some global attributes of the dataset
3169          nc_stat = NF90_PUT_ATT( id_set_prt, NF90_GLOBAL, 'title', &
3170                                  TRIM( run_description_header ) )
3171          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 310 )
3172
3173!
3174!--       Define time coordinate for particles (unlimited dimension)
3175          nc_stat = NF90_DEF_DIM( id_set_prt, 'time', NF90_UNLIMITED, &
3176                                  id_dim_time_prt )
3177          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 311 )
3178
3179          nc_stat = NF90_DEF_VAR( id_set_prt, 'time', NF90_DOUBLE, &
3180                                  id_dim_time_prt, id_var_time_prt )
3181          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 312 )
3182
3183          nc_stat = NF90_PUT_ATT( id_set_prt, id_var_time_prt, 'units', &
3184                                  'seconds' )
3185          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 313 )
3186
3187!
3188!--       Define particle coordinate (maximum particle number)
3189          nc_stat = NF90_DEF_DIM( id_set_prt, 'particle_number', &
3190                                  maximum_number_of_particles, id_dim_prtnum )
3191          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 314 )
3192
3193          nc_stat = NF90_DEF_VAR( id_set_prt, 'particle_number', NF90_DOUBLE, &
3194                                  id_dim_prtnum, id_var_prtnum )
3195          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 315 )
3196
3197          nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prtnum, 'units', &
3198                                  'particle number' )
3199          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 316 )
3200
3201!
3202!--       Define variable which contains the real number of particles in use
3203          nc_stat = NF90_DEF_VAR( id_set_prt, 'real_num_of_prt', NF90_INT, &
3204                                  id_dim_time_prt, id_var_rnop_prt )
3205          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 317 )
3206
3207          nc_stat = NF90_PUT_ATT( id_set_prt, id_var_rnop_prt, 'units', &
3208                                  'particle number' )
3209          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 318 )
3210
3211!
3212!--       Define the variables
3213          DO  i = 1, 17
3214
3215             nc_stat = NF90_DEF_VAR( id_set_prt, prt_var_names(i),         &
3216                                     nc_precision(8),                      &
3217                                     (/ id_dim_prtnum, id_dim_time_prt /), &
3218                                     id_var_prt(i) )
3219             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 319 )
3220
3221             nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prt(i), &
3222                                     'long_name', TRIM( prt_var_names(i) ) )
3223             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 320 )
3224
3225             nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prt(i), &
3226                                     'units', TRIM( prt_var_units(i) ) )
3227             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 321 )
3228
3229          ENDDO
3230
3231!
3232!--       Leave NetCDF define mode
3233          nc_stat = NF90_ENDDEF( id_set_prt )
3234          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 322 )
3235
3236
3237       CASE ( 'pt_ext' )
3238
3239!
3240!--       Get the id of the time coordinate (unlimited coordinate) and its
3241!--       last index on the file. The next time level is prt..count+1.
3242!--       The current time must be larger than the last output time
3243!--       on the file.
3244          nc_stat = NF90_INQ_VARID( id_set_prt, 'time', id_var_time_prt )
3245          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 323 )
3246
3247          nc_stat = NF90_INQUIRE_VARIABLE( id_set_prt, id_var_time_prt, &
3248                                           dimids = id_dim_time_old )
3249          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 324 )
3250          id_dim_time_prt = id_dim_time_old(1)
3251
3252          nc_stat = NF90_INQUIRE_DIMENSION( id_set_prt, id_dim_time_prt, &
3253                                            len = prt_time_count )
3254          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 325 )
3255
3256          nc_stat = NF90_GET_VAR( id_set_prt, id_var_time_prt,  &
3257                                  last_time_coordinate,         &
3258                                  start = (/ prt_time_count /), &
3259                                  count = (/ 1 /) )
3260          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 326 )
3261
3262          IF ( last_time_coordinate(1) >= simulated_time )  THEN
3263             PRINT*, '+++ WARNING: NetCDF file for particles from previous ' //&
3264                                   'run found,'
3265             PRINT*, '             but this file cannot be extended becaus' // &
3266                                   'e the current output time'
3267             PRINT*, '             is less or equal than the last output t' // &
3268                                   'ime on this file.'
3269             PRINT*, '             New file is created instead.'
3270             prt_time_count = 0
3271             extend = .FALSE.
3272             RETURN
3273          ENDIF
3274
3275!
3276!--       Dataset seems to be extendable.
3277!--       Now get the variable ids.
3278          nc_stat = NF90_INQ_VARID( id_set_prt, 'real_num_of_prt', &
3279                                    id_var_rnop_prt )
3280          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 327 )
3281
3282          DO  i = 1, 17
3283
3284             nc_stat = NF90_INQ_VARID( id_set_prt, prt_var_names(i), &
3285                                       id_var_prt(i) )
3286             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 328 )
3287
3288          ENDDO
3289
3290          IF ( myid == 0 )  THEN
3291             PRINT*, '*** NetCDF file for particles from previous run found.'
3292             PRINT*, '    This file will be extended.'
3293          ENDIF
3294
3295
3296       CASE ( 'ps_new' )
3297
3298!
3299!--       Define some global attributes of the dataset
3300          nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
3301                                  TRIM( run_description_header ) )
3302          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 396 )
3303
3304!
3305!--       Define time coordinate for particle time series (unlimited dimension)
3306          nc_stat = NF90_DEF_DIM( id_set_pts, 'time', NF90_UNLIMITED, &
3307                                  id_dim_time_pts )
3308          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 397 )
3309
3310          nc_stat = NF90_DEF_VAR( id_set_pts, 'time', NF90_DOUBLE, &
3311                                  id_dim_time_pts, id_var_time_pts )
3312          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 398 )
3313
3314          nc_stat = NF90_PUT_ATT( id_set_pts, id_var_time_pts, 'units', &
3315                                  'seconds')
3316          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 399 )
3317
3318!
3319!--       Define the variables. If more than one particle group is defined,
3320!--       define seperate variables for each group
3321          var_list = ';'
3322          DO  i = 1, dopts_num
3323
3324!
3325!--          First, remove those characters not allowed by NetCDF
3326             netcdf_var_name = dopts_label(i)
3327             CALL clean_netcdf_varname( netcdf_var_name )
3328
3329             DO  j = 0, number_of_particle_groups
3330
3331                IF ( j == 0 )  THEN
3332                   suffix1 = ''
3333                ELSE
3334                   WRITE ( suffix1, '(''_'',I2.2)' )  j
3335                ENDIF
3336
3337                nc_stat = NF90_DEF_VAR( id_set_pts,                         &
3338                                        TRIM( netcdf_var_name ) // suffix1, &
3339                                        nc_precision(6),                    &
3340                                        (/ id_dim_time_pts /),              &
3341                                        id_var_dopts(i,j) )
3342                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 400 )
3343
3344                IF ( j == 0 )  THEN
3345                   nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
3346                                           'long_name', &
3347                                           TRIM( dopts_label(i) ) )
3348                ELSE
3349                   nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
3350                                           'long_name', &
3351                                           TRIM( dopts_label(i) ) // ' PG ' // &
3352                                           suffix1(2:3) )
3353                ENDIF
3354                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 401 )
3355
3356                nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
3357                                        'units', TRIM( dopts_unit(i) ) )
3358                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 402 )
3359
3360                var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
3361                           suffix1 // ';'
3362
3363                IF ( number_of_particle_groups == 1 )  EXIT
3364
3365             ENDDO
3366
3367          ENDDO
3368
3369!
3370!--       Write the list of variables as global attribute (this is used by
3371!--       restart runs)
3372          nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
3373                                  var_list )
3374          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 403 )
3375
3376!
3377!--       Leave NetCDF define mode
3378          nc_stat = NF90_ENDDEF( id_set_pts )
3379          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 404 )
3380
3381
3382       CASE ( 'ps_ext' )
3383
3384!
3385!--       Get the list of variables and compare with the actual run.
3386!--       First var_list_old has to be reset, since GET_ATT does not assign
3387!--       trailing blanks.
3388          var_list_old = ' '
3389          nc_stat = NF90_GET_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
3390                                  var_list_old )
3391          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 405 )
3392
3393          var_list = ';'
3394          i = 1
3395          DO  i = 1, dopts_num
3396
3397             netcdf_var_name = dopts_label(i)
3398             CALL clean_netcdf_varname( netcdf_var_name )
3399
3400             DO  j = 0, number_of_particle_groups
3401
3402                IF ( j == 0 )  THEN
3403                   suffix1 = ''
3404                ELSE
3405                   WRITE ( suffix1, '(''_'',I2.2)' )  j
3406                ENDIF
3407
3408                var_list = TRIM(var_list) // TRIM(netcdf_var_name) // &
3409                           suffix1 // ';'
3410
3411                IF ( number_of_particle_groups == 1 )  EXIT
3412
3413             ENDDO
3414
3415          ENDDO
3416
3417          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
3418             PRINT*, '+++ WARNING: NetCDF file for particle time series ' //&
3419                                   'from previuos run found,'
3420             PRINT*, '             but this file cannot be extended due to' // &
3421                                   ' variable mismatch.'
3422             PRINT*, '             New file is created instead.'
3423             extend = .FALSE.
3424             RETURN
3425          ENDIF
3426
3427!
3428!--       Get the id of the time coordinate (unlimited coordinate) and its
3429!--       last index on the file. The next time level is dots..count+1.
3430!--       The current time must be larger than the last output time
3431!--       on the file.
3432          nc_stat = NF90_INQ_VARID( id_set_pts, 'time', id_var_time_pts )
3433          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 406 )
3434
3435          nc_stat = NF90_INQUIRE_VARIABLE( id_set_pts, id_var_time_pts, &
3436                                           dimids = id_dim_time_old )
3437          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 407 )
3438          id_dim_time_pts = id_dim_time_old(1)
3439
3440          nc_stat = NF90_INQUIRE_DIMENSION( id_set_pts, id_dim_time_pts, &
3441                                            len = dopts_time_count )
3442          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 408 )
3443
3444          nc_stat = NF90_GET_VAR( id_set_pts, id_var_time_pts,    &
3445                                  last_time_coordinate,           &
3446                                  start = (/ dopts_time_count /), &
3447                                  count = (/ 1 /) )
3448          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 409 )
3449
3450          IF ( last_time_coordinate(1) >= simulated_time )  THEN
3451             PRINT*, '+++ WARNING: NetCDF file for time series from' // &
3452                                   ' previuos run found,'
3453             PRINT*, '             but this file cannot be extended becaus' // &
3454                                   'e the current output time'
3455             PRINT*, '             is less or equal than the last output t' // &
3456                                   'ime on this file.'
3457             PRINT*, '             New file is created instead.'
3458             dopts_time_count = 0
3459             extend = .FALSE.
3460             RETURN
3461          ENDIF
3462
3463!
3464!--       Dataset seems to be extendable.
3465!--       Now get the variable ids
3466          i = 1
3467          DO  i = 1, dopts_num
3468 
3469             netcdf_var_name_base = dopts_label(i)
3470             CALL clean_netcdf_varname( netcdf_var_name_base )
3471
3472             DO  j = 0, number_of_particle_groups
3473
3474                IF ( j == 0 )  THEN
3475                   suffix1 = ''
3476                ELSE
3477                   WRITE ( suffix1, '(''_'',I2.2)' )  j
3478                ENDIF
3479
3480                netcdf_var_name = TRIM( netcdf_var_name_base ) // suffix1
3481
3482                nc_stat = NF90_INQ_VARID( id_set_pts, netcdf_var_name, &
3483                                          id_var_dopts(i,j) )
3484                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 410 )
3485
3486                IF ( number_of_particle_groups == 1 )  EXIT
3487
3488             ENDDO
3489
3490          ENDDO
3491
3492!
3493!--       Change the title attribute on file
3494          nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
3495                                  TRIM( run_description_header ) )
3496          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 411 )
3497          PRINT*, '*** NetCDF file for particle time series from previous ', &
3498                       'run found.'
3499          PRINT*, '    This file will be extended.'
3500
3501
3502       CASE DEFAULT
3503
3504          PRINT*, '+++ define_netcdf_header: mode "', mode, '" not supported'
3505
3506    END SELECT
3507
3508#endif
3509 END SUBROUTINE define_netcdf_header
3510
3511
3512
3513 SUBROUTINE handle_netcdf_error( errno )
3514#if defined( __netcdf )
3515
3516!------------------------------------------------------------------------------!
3517!
3518! Description:
3519! ------------
3520! Prints out a text message corresponding to the current status.
3521!------------------------------------------------------------------------------!
3522
3523    USE netcdf
3524    USE netcdf_control
3525    USE pegrid
3526
3527    IMPLICIT NONE
3528
3529    INTEGER ::  errno
3530
3531    IF ( nc_stat /= NF90_NOERR )  THEN
3532       PRINT*, '+++ netcdf error ', errno,': ', TRIM( NF90_STRERROR( nc_stat ) )
3533#if defined( __parallel )
3534       CALL MPI_ABORT( comm2d, 9999, ierr )
3535#else
3536       CALL local_stop
3537#endif
3538    ENDIF
3539
3540#endif
3541 END SUBROUTINE handle_netcdf_error
3542
3543
3544
3545 SUBROUTINE clean_netcdf_varname( string )
3546#if defined( __netcdf )
3547
3548!------------------------------------------------------------------------------!
3549!
3550! Description:
3551! ------------
3552! Replace those characters in string which are not allowed by NetCDF.
3553!------------------------------------------------------------------------------!
3554
3555    USE netcdf_control
3556
3557    IMPLICIT NONE
3558
3559    CHARACTER (LEN=10), INTENT(INOUT) ::  string
3560
3561    INTEGER ::  i, ic
3562
3563    DO  i = 1, replace_num
3564       DO
3565          ic = INDEX( string, replace_char(i) )
3566          IF ( ic /= 0 )  THEN
3567             string(ic:ic) = replace_by(i)
3568          ELSE
3569             EXIT
3570          ENDIF
3571       ENDDO
3572    ENDDO
3573
3574#endif
3575 END SUBROUTINE clean_netcdf_varname
Note: See TracBrowser for help on using the repository browser.