Home

Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |