Home

Context Navigation

netcdf.f90 @ 403

Last change on this file since 403 was 392, checked in by raasch, 15 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 franke $
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
3725	SUBROUTINE 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.

Context Navigation

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

Download in other formats: