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

Last change on this file since 392 was 392, checked in by raasch, 12 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:


lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:


Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

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