Home

Context Navigation

source: palm/trunk/SOURCE/check_parameters.f90 @ 358

Last change on this file since 358 was 354, checked in by maronga, 15 years ago
xy cross section output of sensible and latent heatflux is now available
Property svn:keywords set to `Id`
File size: 117.9 KB

Line
1	SUBROUTINE check_parameters
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	! initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill'
7	! Check for illegal entries in section_xy\|xz\|yz that exceed nz+1\|ny+1\|nx+1
8	! Coupling with independent precursor runs.
9	! Check particle_color, particle_dvrpsize, color_interval, dvrpsize_interval
10	! Bugfix: pressure included for profile output
11	! Check pressure gradient conditions
12	! topography_grid_convention moved from user_check_parameters
13	! 'single_street_canyon'
14	!
15	! Former revisions:
16	! -----------------
17	! $Id: check_parameters.f90 354 2009-07-13 13:14:25Z heinze $
18	!
19	! 222 2009-01-12 16:04:16Z letzel
20	! +user_check_parameters
21	! Output of messages replaced by message handling routine.
22	! Implementation of an MPI-1 coupling: replaced myid with target_id,
23	! deleted __mpi2 directives
24	! Check that PALM is called with mrun -K parallel for coupling
25	!
26	! 197 2008-09-16 15:29:03Z raasch
27	! Bug fix: Construction of vertical profiles when 10 gradients have been
28	! specified in the parameter list (ug, vg, pt, q, sa, lad)
29	!
30	! Strict grid matching along z is not needed for mg-solver.
31	! Leaf area density (LAD) explicitly set to its surface value at k=0
32	! Case of reading data for recycling included in initializing_actions,
33	! check of turbulent_inflow and calculation of recycling_plane.
34	! q*2 profile added
35	!
36	! 138 2007-11-28 10:03:58Z letzel
37	! Plant canopy added
38	! Allow new case bc_uv_t = 'dirichlet_0' for channel flow.
39	! Multigrid solver allows topography, checking of dt_sort_particles
40	! Bugfix: initializing u_init and v_init in case of ocean runs
41	!
42	! 109 2007-08-28 15:26:47Z letzel
43	! Check coupling_mode and set default (obligatory) values (like boundary
44	! conditions for temperature and fluxes) in case of coupled runs.
45	! +profiles for wp and w"e
46	! Bugfix: Error message concerning output of particle concentration (pc)
47	! modified
48	! More checks and more default values for coupled runs
49	! allow data_output_pr= q, wq, w"q", wq for humidity = .T. (instead of
50	! cloud_physics = .T.)
51	! Rayleigh damping for ocean fixed.
52	! Check and, if necessary, set default value for dt_coupling
53	!
54	! 97 2007-06-21 08:23:15Z raasch
55	! Initial salinity profile is calculated, salinity boundary conditions are
56	! checked,
57	! z_max_do1d is checked only in case of ocean = .f.,
58	! +initial temperature and geostrophic velocity profiles for the ocean version,
59	! use_pt_reference renamed use_reference
60	!
61	! 89 2007-05-25 12:08:31Z raasch
62	! Check for user-defined profiles
63	!
64	! 75 2007-03-22 09:54:05Z raasch
65	! "by_user" allowed as initializing action, -data_output_ts,
66	! leapfrog with non-flat topography not allowed any more, loop_optimization
67	! and pt_reference are checked, moisture renamed humidity,
68	! output of precipitation amount/rate and roughnes length + check
69	! possible negative humidities are avoided in initial profile,
70	! dirichlet/neumann changed to dirichlet/radiation, etc.,
71	! revision added to run_description_header
72	!
73	! 20 2007-02-26 00:12:32Z raasch
74	! Temperature and humidity gradients at top are now calculated for nzt+1,
75	! top_heatflux and respective boundary condition bc_pt_t is checked
76	!
77	! RCS Log replace by Id keyword, revision history cleaned up
78	!
79	! Revision 1.61 2006/08/04 14:20:25 raasch
80	! do2d_unit and do3d_unit now defined as 2d-arrays, check of
81	! use_upstream_for_tke, default value for dt_dopts,
82	! generation of file header moved from routines palm and header to here
83	!
84	! Revision 1.1 1997/08/26 06:29:23 raasch
85	! Initial revision
86	!
87	!
88	! Description:
89	! ------------
90	! Check control parameters and deduce further quantities.
91	!------------------------------------------------------------------------------!
92
93	USE arrays_3d
94	USE constants
95	USE control_parameters
96	USE dvrp_variables
97	USE grid_variables
98	USE indices
99	USE model_1d
100	USE netcdf_control
101	USE particle_attributes
102	USE pegrid
103	USE profil_parameter
104	USE statistics
105	USE transpose_indices
106
107	IMPLICIT NONE
108
109	CHARACTER (LEN=1) :: sq
110	CHARACTER (LEN=6) :: var
111	CHARACTER (LEN=7) :: unit
112	CHARACTER (LEN=8) :: date
113	CHARACTER (LEN=10) :: time
114	CHARACTER (LEN=40) :: coupling_string
115	CHARACTER (LEN=100) :: action
116
117	INTEGER :: i, ilen, intervals, iremote = 0, iter, j, k, nnxh, nnyh, &
118	position, prec
119	LOGICAL :: found, ldum
120	REAL :: gradient, maxn, maxp, remote = 0.0, &
121	simulation_time_since_reference
122
123	!
124	!-- Warning, if host is not set
125	IF ( host(1:1) == ' ' ) THEN
126	message_string = '"host" is not set. Please check that environment ' // &
127	'variable "localhost" & is set before running PALM'
128	CALL message( 'check_parameters', 'PA0001', 0, 0, 0, 6, 0 )
129	ENDIF
130
131	!
132	!-- Check the coupling mode
133	IF ( coupling_mode /= 'uncoupled' .AND. &
134	coupling_mode /= 'atmosphere_to_ocean' .AND. &
135	coupling_mode /= 'ocean_to_atmosphere' ) THEN
136	message_string = 'illegal coupling mode: ' // TRIM( coupling_mode )
137	CALL message( 'check_parameters', 'PA0002', 1, 2, 0, 6, 0 )
138	ENDIF
139
140	!
141	!-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
142	IF ( coupling_mode /= 'uncoupled' ) THEN
143
144	IF ( dt_coupling == 9999999.9 ) THEN
145	message_string = 'dt_coupling is not set but required for coup' // &
146	'ling mode "' // TRIM( coupling_mode ) // '"'
147	CALL message( 'check_parameters', 'PA0003', 1, 2, 0, 6, 0 )
148	ENDIF
149
150	#if defined( __parallel )
151	CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, &
152	ierr )
153	CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, &
154	status, ierr )
155	IF ( dt_coupling /= remote ) THEN
156	WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
157	'": dt_coupling = ', dt_coupling, '& is not equal to ', &
158	'dt_coupling_remote = ', remote
159	CALL message( 'check_parameters', 'PA0004', 1, 2, 0, 6, 0 )
160	ENDIF
161	IF ( dt_coupling <= 0.0 ) THEN
162	CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
163	CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, &
164	status, ierr )
165	dt_coupling = MAX( dt_max, remote )
166	WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
167	'": dt_coupling <= 0.0 & is not allowed and is reset to ', &
168	'MAX(dt_max(A,O)) = ', dt_coupling
169	CALL message( 'check_parameters', 'PA0005', 0, 1, 0, 6, 0 )
170	ENDIF
171
172	CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
173	ierr )
174	CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, &
175	status, ierr )
176	IF ( restart_time /= remote ) THEN
177	WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
178	'": restart_time = ', restart_time, '& is not equal to ', &
179	'restart_time_remote = ', remote
180	CALL message( 'check_parameters', 'PA0006', 1, 2, 0, 6, 0 )
181	ENDIF
182
183	CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, &
184	ierr )
185	CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, &
186	status, ierr )
187	IF ( dt_restart /= remote ) THEN
188	WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
189	'": dt_restart = ', dt_restart, '& is not equal to ', &
190	'dt_restart_remote = ', remote
191	CALL message( 'check_parameters', 'PA0007', 1, 2, 0, 6, 0 )
192	ENDIF
193
194	simulation_time_since_reference = end_time - coupling_start_time
195	CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &
196	14, comm_inter, ierr )
197	CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, &
198	status, ierr )
199	IF ( simulation_time_since_reference /= remote ) THEN
200	WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
201	'": simulation_time_since_reference = ', &
202	simulation_time_since_reference, '& is not equal to ', &
203	'simulation_time_since_reference_remote = ', remote
204	CALL message( 'check_parameters', 'PA0008', 1, 2, 0, 6, 0 )
205	ENDIF
206
207	CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
208	CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, &
209	status, ierr )
210	IF ( dx /= remote ) THEN
211	WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
212	'": dx = ', dx, '& is not equal to dx_remote = ', remote
213	CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
214	ENDIF
215
216	CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
217	CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, &
218	status, ierr )
219	IF ( dy /= remote ) THEN
220	WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
221	'": dy = ', dy, '& is not equal to dy_remote = ', remote
222	CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
223	ENDIF
224
225	CALL MPI_SEND( nx, 1, MPI_INTEGER, target_id, 17, comm_inter, ierr )
226	CALL MPI_RECV( iremote, 1, MPI_INTEGER, target_id, 17, comm_inter, &
227	status, ierr )
228	IF ( nx /= iremote ) THEN
229	WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
230	'": nx = ', nx, '& is not equal to nx_remote = ', iremote
231	CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
232	ENDIF
233
234	CALL MPI_SEND( ny, 1, MPI_INTEGER, target_id, 18, comm_inter, ierr )
235	CALL MPI_RECV( iremote, 1, MPI_INTEGER, target_id, 18, comm_inter, &
236	status, ierr )
237	IF ( ny /= iremote ) THEN
238	WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
239	'": ny = ', ny, '& is not equal to ny_remote = ', iremote
240	CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
241	ENDIF
242	#else
243	WRITE( message_string, * ) 'coupling requires PALM to be called with', &
244	' ''mrun -K parallel'''
245	CALL message( 'check_parameters', 'PA0141', 1, 2, 0, 6, 0 )
246	#endif
247	ENDIF
248
249	#if defined( __parallel )
250	!
251	!-- Exchange via intercommunicator
252	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
253	CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter, &
254	ierr )
255	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
256	CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19, &
257	comm_inter, status, ierr )
258	ENDIF
259	#endif
260
261
262	!
263	!-- Generate the file header which is used as a header for most of PALM's
264	!-- output files
265	CALL DATE_AND_TIME( date, time )
266	run_date = date(7:8)//'-'//date(5:6)//'-'//date(3:4)
267	run_time = time(1:2)//':'//time(3:4)//':'//time(5:6)
268	IF ( coupling_mode == 'uncoupled' ) THEN
269	coupling_string = ''
270	ELSEIF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
271	coupling_string = ' coupled (atmosphere)'
272	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
273	coupling_string = ' coupled (ocean)'
274	ENDIF
275
276	WRITE ( run_description_header, &
277	'(A,2X,A,2X,A,A,A,I2.2,A,2X,A,A,2X,A,1X,A)' ) &
278	TRIM( version ), TRIM( revision ), 'run: ', &
279	TRIM( run_identifier ), '.', runnr, TRIM( coupling_string ), &
280	'host: ', TRIM( host ), run_date, run_time
281
282	!
283	!-- Check the general loop optimization method
284	IF ( loop_optimization == 'default' ) THEN
285	IF ( host(1:3) == 'nec' ) THEN
286	loop_optimization = 'vector'
287	ELSE
288	loop_optimization = 'cache'
289	ENDIF
290	ENDIF
291	IF ( loop_optimization /= 'noopt' .AND. loop_optimization /= 'cache' &
292	.AND. loop_optimization /= 'vector' ) THEN
293	message_string = 'illegal value given for loop_optimization: "' // &
294	TRIM( loop_optimization ) // '"'
295	CALL message( 'check_parameters', 'PA0013', 1, 2, 0, 6, 0 )
296	ENDIF
297
298	!
299	!-- Check topography setting (check for illegal parameter combinations)
300	IF ( topography /= 'flat' ) THEN
301	action = ' '
302	IF ( scalar_advec /= 'pw-scheme' ) THEN
303	WRITE( action, '(A,A)' ) 'scalar_advec = ', scalar_advec
304	ENDIF
305	IF ( momentum_advec /= 'pw-scheme' ) THEN
306	WRITE( action, '(A,A)' ) 'momentum_advec = ', momentum_advec
307	ENDIF
308	IF ( timestep_scheme(1:8) == 'leapfrog' ) THEN
309	WRITE( action, '(A,A)' ) 'timestep_scheme = ', timestep_scheme
310	ENDIF
311	IF ( psolver == 'sor' ) THEN
312	WRITE( action, '(A,A)' ) 'psolver = ', psolver
313	ENDIF
314	IF ( sloping_surface ) THEN
315	WRITE( action, '(A)' ) 'sloping surface = .TRUE.'
316	ENDIF
317	IF ( galilei_transformation ) THEN
318	WRITE( action, '(A)' ) 'galilei_transformation = .TRUE.'
319	ENDIF
320	IF ( cloud_physics ) THEN
321	WRITE( action, '(A)' ) 'cloud_physics = .TRUE.'
322	ENDIF
323	IF ( cloud_droplets ) THEN
324	WRITE( action, '(A)' ) 'cloud_droplets = .TRUE.'
325	ENDIF
326	IF ( humidity ) THEN
327	WRITE( action, '(A)' ) 'humidity = .TRUE.'
328	ENDIF
329	IF ( .NOT. prandtl_layer ) THEN
330	WRITE( action, '(A)' ) 'prandtl_layer = .FALSE.'
331	ENDIF
332	IF ( action /= ' ' ) THEN
333	message_string = 'a non-flat topography does not allow ' // &
334	TRIM( action )
335	CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
336	ENDIF
337	!
338	!-- In case of non-flat topography, check whether the convention how to
339	!-- define the topography grid has been set correctly, or whether the default
340	!-- is applicable. If this is not possible, abort.
341	IF ( TRIM( topography_grid_convention ) == ' ' ) THEN
342	IF ( TRIM( topography ) /= 'single_building' .AND. &
343	TRIM( topography ) /= 'single_street_canyon' .AND. &
344	TRIM( topography ) /= 'read_from_file' ) THEN
345	!-- The default value is not applicable here, because it is only valid
346	!-- for the two standard cases 'single_building' and 'read_from_file'
347	!-- defined in init_grid.
348	WRITE( message_string, * ) &
349	'The value for "topography_grid_convention" ', &
350	'is not set. Its default value is & only valid for ', &
351	'"topography" = ''single_building'', ', &
352	'''single_street_canyon'' & or ''read_from_file''.', &
353	' & Choose ''cell_edge'' or ''cell_center''.'
354	CALL message( 'user_check_parameters', 'PA0239', 1, 2, 0, 6, 0 )
355	ELSE
356	!-- The default value is applicable here.
357	!-- Set convention according to topography.
358	IF ( TRIM( topography ) == 'single_building' .OR. &
359	TRIM( topography ) == 'single_street_canyon' ) THEN
360	topography_grid_convention = 'cell_edge'
361	ELSEIF ( TRIM( topography ) == 'read_from_file' ) THEN
362	topography_grid_convention = 'cell_center'
363	ENDIF
364	ENDIF
365	ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND. &
366	TRIM( topography_grid_convention ) /= 'cell_center' ) THEN
367	WRITE( message_string, * ) &
368	'The value for "topography_grid_convention" is ', &
369	'not recognized. & Choose ''cell_edge'' or ''cell_center''.'
370	CALL message( 'user_check_parameters', 'PA0240', 1, 2, 0, 6, 0 )
371	ENDIF
372
373	ENDIF
374
375	!
376	!-- Check ocean setting
377	IF ( ocean ) THEN
378
379	action = ' '
380	IF ( timestep_scheme(1:8) == 'leapfrog' ) THEN
381	WRITE( action, '(A,A)' ) 'timestep_scheme = ', timestep_scheme
382	ENDIF
383	IF ( momentum_advec == 'ups-scheme' ) THEN
384	WRITE( action, '(A,A)' ) 'momentum_advec = ', momentum_advec
385	ENDIF
386	IF ( action /= ' ' ) THEN
387	message_string = 'ocean = .T. does not allow ' // TRIM( action )
388	CALL message( 'check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
389	ENDIF
390
391	ELSEIF ( TRIM( coupling_mode ) == 'uncoupled' .AND. &
392	TRIM( coupling_char ) == '_O' ) THEN
393
394	!
395	!-- Check whether an (uncoupled) atmospheric run has been declared as an
396	!-- ocean run (this setting is done via mrun-option -y)
397
398	message_string = 'ocean = .F. does not allow coupling_char = "' // &
399	TRIM( coupling_char ) // '" set by mrun-option "-y"'
400	CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
401
402	ENDIF
403
404	!
405	!-- Check whether there are any illegal values
406	!-- Pressure solver:
407	IF ( psolver /= 'poisfft' .AND. psolver /= 'poisfft_hybrid' .AND. &
408	psolver /= 'sor' .AND. psolver /= 'multigrid' ) THEN
409	message_string = 'unknown solver for perturbation pressure: psolver' // &
410	' = "' // TRIM( psolver ) // '"'
411	CALL message( 'check_parameters', 'PA0016', 1, 2, 0, 6, 0 )
412	ENDIF
413
414	#if defined( __parallel )
415	IF ( psolver == 'poisfft_hybrid' .AND. pdims(2) /= 1 ) THEN
416	message_string = 'psolver = "' // TRIM( psolver ) // '" only works ' // &
417	'for a 1d domain-decomposition along x & please do' // &
418	' not set npey/=1 in the parameter file'
419	CALL message( 'check_parameters', 'PA0017', 1, 2, 0, 6, 0 )
420	ENDIF
421	IF ( psolver == 'poisfft_hybrid' .AND. &
422	( nxra > nxr .OR. nyna > nyn .OR. nza > nz ) .OR. &
423	psolver == 'multigrid' .AND. &
424	( nxra > nxr .OR. nyna > nyn ) ) THEN
425	message_string = 'psolver = "' // TRIM( psolver ) // '" does not ' // &
426	'work for subdomains with unequal size & please ' // &
427	'set grid_matching = ''strict'' in the parameter file'
428	CALL message( 'check_parameters', 'PA0018', 1, 2, 0, 6, 0 )
429	ENDIF
430	#else
431	IF ( psolver == 'poisfft_hybrid' ) THEN
432	message_string = 'psolver = "' // TRIM( psolver ) // '" only works' // &
433	' for a parallel environment'
434	CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
435	ENDIF
436	#endif
437
438	IF ( psolver == 'multigrid' ) THEN
439	IF ( cycle_mg == 'w' ) THEN
440	gamma_mg = 2
441	ELSEIF ( cycle_mg == 'v' ) THEN
442	gamma_mg = 1
443	ELSE
444	message_string = 'unknown multigrid cycle: cycle_mg = "' // &
445	TRIM( cycle_mg ) // '"'
446	CALL message( 'check_parameters', 'PA0020', 1, 2, 0, 6, 0 )
447	ENDIF
448	ENDIF
449
450	IF ( fft_method /= 'singleton-algorithm' .AND. &
451	fft_method /= 'temperton-algorithm' .AND. &
452	fft_method /= 'system-specific' ) THEN
453	message_string = 'unknown fft-algorithm: fft_method = "' // &
454	TRIM( fft_method ) // '"'
455	CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
456	ENDIF
457
458	!
459	!-- Advection schemes:
460	IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ups-scheme' ) &
461	THEN
462	message_string = 'unknown advection scheme: momentum_advec = "' // &
463	TRIM( momentum_advec ) // '"'
464	CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
465	ENDIF
466	IF ( ( momentum_advec == 'ups-scheme' .OR. scalar_advec == 'ups-scheme' )&
467	.AND. timestep_scheme /= 'euler' ) THEN
468	message_string = 'momentum_advec = "' // TRIM( momentum_advec ) // &
469	'" is not allowed with timestep_scheme = "' // &
470	TRIM( timestep_scheme ) // '"'
471	CALL message( 'check_parameters', 'PA0023', 1, 2, 0, 6, 0 )
472	ENDIF
473
474	IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'bc-scheme' .AND.&
475	scalar_advec /= 'ups-scheme' ) THEN
476	message_string = 'unknown advection scheme: scalar_advec = "' // &
477	TRIM( scalar_advec ) // '"'
478	CALL message( 'check_parameters', 'PA0024', 1, 2, 0, 6, 0 )
479	ENDIF
480
481	IF ( use_sgs_for_particles .AND. .NOT. use_upstream_for_tke ) THEN
482	use_upstream_for_tke = .TRUE.
483	message_string = 'use_upstream_for_tke set .TRUE. because ' // &
484	'use_sgs_for_particles = .TRUE.'
485	CALL message( 'check_parameters', 'PA0025', 0, 1, 0, 6, 0 )
486	ENDIF
487
488	IF ( use_upstream_for_tke .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
489	message_string = 'use_upstream_for_tke = .TRUE. not allowed with ' // &
490	'timestep_scheme = "' // TRIM( timestep_scheme ) // '"'
491	CALL message( 'check_parameters', 'PA0026', 1, 2, 0, 6, 0 )
492	ENDIF
493
494	!
495	!-- Timestep schemes:
496	SELECT CASE ( TRIM( timestep_scheme ) )
497
498	CASE ( 'euler' )
499	intermediate_timestep_count_max = 1
500	asselin_filter_factor = 0.0
501
502	CASE ( 'leapfrog', 'leapfrog+euler' )
503	intermediate_timestep_count_max = 1
504
505	CASE ( 'runge-kutta-2' )
506	intermediate_timestep_count_max = 2
507	asselin_filter_factor = 0.0
508
509	CASE ( 'runge-kutta-3' )
510	intermediate_timestep_count_max = 3
511	asselin_filter_factor = 0.0
512
513	CASE DEFAULT
514	message_string = 'unknown timestep scheme: timestep_scheme = "' // &
515	TRIM( timestep_scheme ) // '"'
516	CALL message( 'check_parameters', 'PA0027', 1, 2, 0, 6, 0 )
517
518	END SELECT
519
520	IF ( scalar_advec == 'ups-scheme' .AND. timestep_scheme(1:5) == 'runge' )&
521	THEN
522	message_string = 'scalar advection scheme "' // TRIM( scalar_advec ) // &
523	'" & does not work with timestep_scheme "' // &
524	TRIM( timestep_scheme ) // '"'
525	CALL message( 'check_parameters', 'PA0028', 1, 2, 0, 6, 0 )
526	ENDIF
527
528	IF ( momentum_advec /= 'pw-scheme' .AND. timestep_scheme(1:5) == 'runge' ) &
529	THEN
530	message_string = 'momentum advection scheme "' // &
531	TRIM( momentum_advec ) // '" & does not work with ' // &
532	'timestep_scheme "' // TRIM( timestep_scheme ) // '"'
533	CALL message( 'check_parameters', 'PA0029', 1, 2, 0, 6, 0 )
534	ENDIF
535
536	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
537	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
538	!
539	!-- No restart run: several initialising actions are possible
540	action = initializing_actions
541	DO WHILE ( TRIM( action ) /= '' )
542	position = INDEX( action, ' ' )
543	SELECT CASE ( action(1:position-1) )
544
545	CASE ( 'set_constant_profiles', 'set_1d-model_profiles', &
546	'by_user', 'initialize_vortex', 'initialize_ptanom' )
547	action = action(position+1:)
548
549	CASE DEFAULT
550	message_string = 'initializing_action = "' // &
551	TRIM( action ) // '" unkown or not allowed'
552	CALL message( 'check_parameters', 'PA0030', 1, 2, 0, 6, 0 )
553
554	END SELECT
555	ENDDO
556	ENDIF
557
558	IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .AND. &
559	INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
560	message_string = 'initializing_actions = "set_constant_profiles"' // &
561	' and "set_1d-model_profiles" are not allowed ' // &
562	'simultaneously'
563	CALL message( 'check_parameters', 'PA0031', 1, 2, 0, 6, 0 )
564	ENDIF
565
566	IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .AND. &
567	INDEX( initializing_actions, 'by_user' ) /= 0 ) THEN
568	message_string = 'initializing_actions = "set_constant_profiles"' // &
569	' and "by_user" are not allowed simultaneously'
570	CALL message( 'check_parameters', 'PA0032', 1, 2, 0, 6, 0 )
571	ENDIF
572
573	IF ( INDEX( initializing_actions, 'by_user' ) /= 0 .AND. &
574	INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
575	message_string = 'initializing_actions = "by_user" and ' // &
576	'"set_1d-model_profiles" are not allowed simultaneously'
577	CALL message( 'check_parameters', 'PA0033', 1, 2, 0, 6, 0 )
578	ENDIF
579
580	IF ( cloud_physics .AND. .NOT. humidity ) THEN
581	WRITE( message_string, * ) 'cloud_physics = ', cloud_physics, ' is ', &
582	'not allowed with humidity = ', humidity
583	CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
584	ENDIF
585
586	IF ( precipitation .AND. .NOT. cloud_physics ) THEN
587	WRITE( message_string, * ) 'precipitation = ', precipitation, ' is ', &
588	'not allowed with cloud_physics = ', cloud_physics
589	CALL message( 'check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
590	ENDIF
591
592	IF ( humidity .AND. sloping_surface ) THEN
593	message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' // &
594	'are not allowed simultaneously'
595	CALL message( 'check_parameters', 'PA0036', 1, 2, 0, 6, 0 )
596	ENDIF
597
598	IF ( humidity .AND. scalar_advec == 'ups-scheme' ) THEN
599	message_string = 'UPS-scheme is not implemented for humidity = .TRUE.'
600	CALL message( 'check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
601	ENDIF
602
603	IF ( passive_scalar .AND. humidity ) THEN
604	message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' // &
605	'is not allowed simultaneously'
606	CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
607	ENDIF
608
609	IF ( passive_scalar .AND. scalar_advec == 'ups-scheme' ) THEN
610	message_string = 'UPS-scheme is not implemented for passive_scalar' // &
611	' = .TRUE.'
612	CALL message( 'check_parameters', 'PA0039', 1, 2, 0, 6, 0 )
613	ENDIF
614
615	IF ( grid_matching /= 'strict' .AND. grid_matching /= 'match' ) THEN
616	message_string = 'illegal value "' // TRIM( grid_matching ) // &
617	'" found for parameter grid_matching'
618	CALL message( 'check_parameters', 'PA0040', 1, 2, 0, 6, 0 )
619	ENDIF
620
621	IF ( plant_canopy .AND. ( drag_coefficient == 0.0 ) ) THEN
622	message_string = 'plant_canopy = .TRUE. requires a non-zero drag ' // &
623	'coefficient & given value is drag_coefficient = 0.0'
624	CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
625	ENDIF
626
627	!
628	!-- In case of no model continuation run, check initialising parameters and
629	!-- deduce further quantities
630	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
631
632	!
633	!-- Initial profiles for 1D and 3D model, respectively
634	u_init = ug_surface
635	v_init = vg_surface
636	pt_init = pt_surface
637	IF ( humidity ) q_init = q_surface
638	IF ( ocean ) sa_init = sa_surface
639	IF ( passive_scalar ) q_init = s_surface
640	IF ( plant_canopy ) lad = 0.0
641
642	!
643	!--
644	!-- If required, compute initial profile of the geostrophic wind
645	!-- (component ug)
646	i = 1
647	gradient = 0.0
648
649	IF ( .NOT. ocean ) THEN
650
651	ug_vertical_gradient_level_ind(1) = 0
652	ug(0) = ug_surface
653	DO k = 1, nzt+1
654	IF ( i < 11 ) THEN
655	IF ( ug_vertical_gradient_level(i) < zu(k) .AND. &
656	ug_vertical_gradient_level(i) >= 0.0 ) THEN
657	gradient = ug_vertical_gradient(i) / 100.0
658	ug_vertical_gradient_level_ind(i) = k - 1
659	i = i + 1
660	ENDIF
661	ENDIF
662	IF ( gradient /= 0.0 ) THEN
663	IF ( k /= 1 ) THEN
664	ug(k) = ug(k-1) + dzu(k) * gradient
665	ELSE
666	ug(k) = ug_surface + 0.5 * dzu(k) * gradient
667	ENDIF
668	ELSE
669	ug(k) = ug(k-1)
670	ENDIF
671	ENDDO
672
673	ELSE
674
675	ug_vertical_gradient_level_ind(1) = nzt+1
676	ug(nzt+1) = ug_surface
677	DO k = nzt, 0, -1
678	IF ( i < 11 ) THEN
679	IF ( ug_vertical_gradient_level(i) > zu(k) .AND. &
680	ug_vertical_gradient_level(i) <= 0.0 ) THEN
681	gradient = ug_vertical_gradient(i) / 100.0
682	ug_vertical_gradient_level_ind(i) = k + 1
683	i = i + 1
684	ENDIF
685	ENDIF
686	IF ( gradient /= 0.0 ) THEN
687	IF ( k /= nzt ) THEN
688	ug(k) = ug(k+1) - dzu(k+1) * gradient
689	ELSE
690	ug(k) = ug_surface - 0.5 * dzu(k+1) * gradient
691	ug(k+1) = ug_surface + 0.5 * dzu(k+1) * gradient
692	ENDIF
693	ELSE
694	ug(k) = ug(k+1)
695	ENDIF
696	ENDDO
697
698	ENDIF
699
700	u_init = ug
701
702	!
703	!-- In case of no given gradients for ug, choose a vanishing gradient
704	IF ( ug_vertical_gradient_level(1) == -9999999.9 ) THEN
705	ug_vertical_gradient_level(1) = 0.0
706	ENDIF
707
708	!
709	!--
710	!-- If required, compute initial profile of the geostrophic wind
711	!-- (component vg)
712	i = 1
713	gradient = 0.0
714
715	IF ( .NOT. ocean ) THEN
716
717	vg_vertical_gradient_level_ind(1) = 0
718	vg(0) = vg_surface
719	DO k = 1, nzt+1
720	IF ( i < 11 ) THEN
721	IF ( vg_vertical_gradient_level(i) < zu(k) .AND. &
722	vg_vertical_gradient_level(i) >= 0.0 ) THEN
723	gradient = vg_vertical_gradient(i) / 100.0
724	vg_vertical_gradient_level_ind(i) = k - 1
725	i = i + 1
726	ENDIF
727	ENDIF
728	IF ( gradient /= 0.0 ) THEN
729	IF ( k /= 1 ) THEN
730	vg(k) = vg(k-1) + dzu(k) * gradient
731	ELSE
732	vg(k) = vg_surface + 0.5 * dzu(k) * gradient
733	ENDIF
734	ELSE
735	vg(k) = vg(k-1)
736	ENDIF
737	ENDDO
738
739	ELSE
740
741	vg_vertical_gradient_level_ind(1) = nzt+1
742	vg(nzt+1) = vg_surface
743	DO k = nzt, 0, -1
744	IF ( i < 11 ) THEN
745	IF ( vg_vertical_gradient_level(i) > zu(k) .AND. &
746	vg_vertical_gradient_level(i) <= 0.0 ) THEN
747	gradient = vg_vertical_gradient(i) / 100.0
748	vg_vertical_gradient_level_ind(i) = k + 1
749	i = i + 1
750	ENDIF
751	ENDIF
752	IF ( gradient /= 0.0 ) THEN
753	IF ( k /= nzt ) THEN
754	vg(k) = vg(k+1) - dzu(k+1) * gradient
755	ELSE
756	vg(k) = vg_surface - 0.5 * dzu(k+1) * gradient
757	vg(k+1) = vg_surface + 0.5 * dzu(k+1) * gradient
758	ENDIF
759	ELSE
760	vg(k) = vg(k+1)
761	ENDIF
762	ENDDO
763
764	ENDIF
765
766	v_init = vg
767
768	!
769	!-- In case of no given gradients for vg, choose a vanishing gradient
770	IF ( vg_vertical_gradient_level(1) == -9999999.9 ) THEN
771	vg_vertical_gradient_level(1) = 0.0
772	ENDIF
773
774	!
775	!-- Compute initial temperature profile using the given temperature gradients
776	i = 1
777	gradient = 0.0
778
779	IF ( .NOT. ocean ) THEN
780
781	pt_vertical_gradient_level_ind(1) = 0
782	DO k = 1, nzt+1
783	IF ( i < 11 ) THEN
784	IF ( pt_vertical_gradient_level(i) < zu(k) .AND. &
785	pt_vertical_gradient_level(i) >= 0.0 ) THEN
786	gradient = pt_vertical_gradient(i) / 100.0
787	pt_vertical_gradient_level_ind(i) = k - 1
788	i = i + 1
789	ENDIF
790	ENDIF
791	IF ( gradient /= 0.0 ) THEN
792	IF ( k /= 1 ) THEN
793	pt_init(k) = pt_init(k-1) + dzu(k) * gradient
794	ELSE
795	pt_init(k) = pt_surface + 0.5 * dzu(k) * gradient
796	ENDIF
797	ELSE
798	pt_init(k) = pt_init(k-1)
799	ENDIF
800	ENDDO
801
802	ELSE
803
804	pt_vertical_gradient_level_ind(1) = nzt+1
805	DO k = nzt, 0, -1
806	IF ( i < 11 ) THEN
807	IF ( pt_vertical_gradient_level(i) > zu(k) .AND. &
808	pt_vertical_gradient_level(i) <= 0.0 ) THEN
809	gradient = pt_vertical_gradient(i) / 100.0
810	pt_vertical_gradient_level_ind(i) = k + 1
811	i = i + 1
812	ENDIF
813	ENDIF
814	IF ( gradient /= 0.0 ) THEN
815	IF ( k /= nzt ) THEN
816	pt_init(k) = pt_init(k+1) - dzu(k+1) * gradient
817	ELSE
818	pt_init(k) = pt_surface - 0.5 * dzu(k+1) * gradient
819	pt_init(k+1) = pt_surface + 0.5 * dzu(k+1) * gradient
820	ENDIF
821	ELSE
822	pt_init(k) = pt_init(k+1)
823	ENDIF
824	ENDDO
825
826	ENDIF
827
828	!
829	!-- In case of no given temperature gradients, choose gradient of neutral
830	!-- stratification
831	IF ( pt_vertical_gradient_level(1) == -9999999.9 ) THEN
832	pt_vertical_gradient_level(1) = 0.0
833	ENDIF
834
835	!
836	!-- Store temperature gradient at the top boundary for possible Neumann
837	!-- boundary condition
838	bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
839
840	!
841	!-- If required, compute initial humidity or scalar profile using the given
842	!-- humidity/scalar gradient. In case of scalar transport, initially store
843	!-- values of the scalar parameters on humidity parameters
844	IF ( passive_scalar ) THEN
845	bc_q_b = bc_s_b
846	bc_q_t = bc_s_t
847	q_surface = s_surface
848	q_surface_initial_change = s_surface_initial_change
849	q_vertical_gradient = s_vertical_gradient
850	q_vertical_gradient_level = s_vertical_gradient_level
851	surface_waterflux = surface_scalarflux
852	ENDIF
853
854	IF ( humidity .OR. passive_scalar ) THEN
855
856	i = 1
857	gradient = 0.0
858	q_vertical_gradient_level_ind(1) = 0
859	DO k = 1, nzt+1
860	IF ( i < 11 ) THEN
861	IF ( q_vertical_gradient_level(i) < zu(k) .AND. &
862	q_vertical_gradient_level(i) >= 0.0 ) THEN
863	gradient = q_vertical_gradient(i) / 100.0
864	q_vertical_gradient_level_ind(i) = k - 1
865	i = i + 1
866	ENDIF
867	ENDIF
868	IF ( gradient /= 0.0 ) THEN
869	IF ( k /= 1 ) THEN
870	q_init(k) = q_init(k-1) + dzu(k) * gradient
871	ELSE
872	q_init(k) = q_init(k-1) + 0.5 * dzu(k) * gradient
873	ENDIF
874	ELSE
875	q_init(k) = q_init(k-1)
876	ENDIF
877	!
878	!-- Avoid negative humidities
879	IF ( q_init(k) < 0.0 ) THEN
880	q_init(k) = 0.0
881	ENDIF
882	ENDDO
883
884	!
885	!-- In case of no given humidity gradients, choose zero gradient
886	!-- conditions
887	IF ( q_vertical_gradient_level(1) == -1.0 ) THEN
888	q_vertical_gradient_level(1) = 0.0
889	ENDIF
890
891	!
892	!-- Store humidity gradient at the top boundary for possile Neumann
893	!-- boundary condition
894	bc_q_t_val = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
895
896	ENDIF
897
898	!
899	!-- If required, compute initial salinity profile using the given salinity
900	!-- gradients
901	IF ( ocean ) THEN
902
903	i = 1
904	gradient = 0.0
905
906	sa_vertical_gradient_level_ind(1) = nzt+1
907	DO k = nzt, 0, -1
908	IF ( i < 11 ) THEN
909	IF ( sa_vertical_gradient_level(i) > zu(k) .AND. &
910	sa_vertical_gradient_level(i) <= 0.0 ) THEN
911	gradient = sa_vertical_gradient(i) / 100.0
912	sa_vertical_gradient_level_ind(i) = k + 1
913	i = i + 1
914	ENDIF
915	ENDIF
916	IF ( gradient /= 0.0 ) THEN
917	IF ( k /= nzt ) THEN
918	sa_init(k) = sa_init(k+1) - dzu(k+1) * gradient
919	ELSE
920	sa_init(k) = sa_surface - 0.5 * dzu(k+1) * gradient
921	sa_init(k+1) = sa_surface + 0.5 * dzu(k+1) * gradient
922	ENDIF
923	ELSE
924	sa_init(k) = sa_init(k+1)
925	ENDIF
926	ENDDO
927
928	ENDIF
929
930	!
931	!-- If required compute the profile of leaf area density used in the plant canopy model
932	IF ( plant_canopy ) THEN
933
934	i = 1
935	gradient = 0.0
936
937	IF ( .NOT. ocean ) THEN
938
939	lad(0) = lad_surface
940
941	lad_vertical_gradient_level_ind(1) = 0
942	DO k = 1, pch_index
943	IF ( i < 11 ) THEN
944	IF ( lad_vertical_gradient_level(i) < zu(k) .AND. &
945	lad_vertical_gradient_level(i) >= 0.0 ) THEN
946	gradient = lad_vertical_gradient(i)
947	lad_vertical_gradient_level_ind(i) = k - 1
948	i = i + 1
949	ENDIF
950	ENDIF
951	IF ( gradient /= 0.0 ) THEN
952	IF ( k /= 1 ) THEN
953	lad(k) = lad(k-1) + dzu(k) * gradient
954	ELSE
955	lad(k) = lad_surface + 0.5 * dzu(k) *gradient
956	ENDIF
957	ELSE
958	lad(k) = lad(k-1)
959	ENDIF
960	ENDDO
961
962	ENDIF
963
964	!
965	!-- In case of no given leaf area density gradients, choose a vanishing gradient
966	IF ( lad_vertical_gradient_level(1) == -9999999.9 ) THEN
967	lad_vertical_gradient_level(1) = 0.0
968	ENDIF
969
970	ENDIF
971
972	ENDIF
973
974	!
975	!-- Compute Coriolis parameter
976	f = 2.0 * omega * SIN( phi / 180.0 * pi )
977	fs = 2.0 * omega * COS( phi / 180.0 * pi )
978
979	!
980	!-- Ocean runs always use reference values in the buoyancy term. Therefore
981	!-- set the reference temperature equal to the surface temperature.
982	IF ( ocean .AND. pt_reference == 9999999.9 ) pt_reference = pt_surface
983
984	!
985	!-- Reference value has to be used in buoyancy terms
986	IF ( pt_reference /= 9999999.9 ) use_reference = .TRUE.
987
988	!
989	!-- Sign of buoyancy/stability terms
990	IF ( ocean ) atmos_ocean_sign = -1.0
991
992	!
993	!-- Ocean version must use flux boundary conditions at the top
994	IF ( ocean .AND. .NOT. use_top_fluxes ) THEN
995	message_string = 'use_top_fluxes must be .TRUE. in ocean version'
996	CALL message( 'check_parameters', 'PA0042', 1, 2, 0, 6, 0 )
997	ENDIF
998
999	!
1000	!-- In case of a given slope, compute the relevant quantities
1001	IF ( alpha_surface /= 0.0 ) THEN
1002	IF ( ABS( alpha_surface ) > 90.0 ) THEN
1003	WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface, &
1004	' ) must be < 90.0'
1005	CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 )
1006	ENDIF
1007	sloping_surface = .TRUE.
1008	cos_alpha_surface = COS( alpha_surface / 180.0 * pi )
1009	sin_alpha_surface = SIN( alpha_surface / 180.0 * pi )
1010	ENDIF
1011
1012	!
1013	!-- Check time step and cfl_factor
1014	IF ( dt /= -1.0 ) THEN
1015	IF ( dt <= 0.0 .AND. dt /= -1.0 ) THEN
1016	WRITE( message_string, * ) 'dt = ', dt , ' <= 0.0'
1017	CALL message( 'check_parameters', 'PA0044', 1, 2, 0, 6, 0 )
1018	ENDIF
1019	dt_3d = dt
1020	dt_fixed = .TRUE.
1021	ENDIF
1022
1023	IF ( cfl_factor <= 0.0 .OR. cfl_factor > 1.0 ) THEN
1024	IF ( cfl_factor == -1.0 ) THEN
1025	IF ( momentum_advec == 'ups-scheme' .OR. &
1026	scalar_advec == 'ups-scheme' ) THEN
1027	cfl_factor = 0.8
1028	ELSE
1029	IF ( timestep_scheme == 'runge-kutta-2' ) THEN
1030	cfl_factor = 0.8
1031	ELSEIF ( timestep_scheme == 'runge-kutta-3' ) THEN
1032	cfl_factor = 0.9
1033	ELSE
1034	cfl_factor = 0.1
1035	ENDIF
1036	ENDIF
1037	ELSE
1038	WRITE( message_string, * ) 'cfl_factor = ', cfl_factor, &
1039	' out of range & 0.0 < cfl_factor <= 1.0 is required'
1040	CALL message( 'check_parameters', 'PA0045', 1, 2, 0, 6, 0 )
1041	ENDIF
1042	ENDIF
1043
1044	!
1045	!-- Store simulated time at begin
1046	simulated_time_at_begin = simulated_time
1047
1048	!
1049	!-- Store reference time for coupled runs and change the coupling flag,
1050	!-- if ...
1051	IF ( simulated_time == 0.0 ) THEN
1052	IF ( coupling_start_time == 0.0 ) THEN
1053	time_since_reference_point = 0.0
1054	ELSEIF ( time_since_reference_point < 0.0 ) THEN
1055	run_coupled = .FALSE.
1056	ENDIF
1057	ENDIF
1058
1059	!
1060	!-- Set wind speed in the Galilei-transformed system
1061	IF ( galilei_transformation ) THEN
1062	IF ( use_ug_for_galilei_tr .AND. &
1063	ug_vertical_gradient_level(1) == 0.0 .AND. &
1064	vg_vertical_gradient_level(1) == 0.0 ) THEN
1065	u_gtrans = ug_surface
1066	v_gtrans = vg_surface
1067	ELSEIF ( use_ug_for_galilei_tr .AND. &
1068	ug_vertical_gradient_level(1) /= 0.0 ) THEN
1069	message_string = 'baroclinicity (ug) not allowed simultaneously' // &
1070	' with galilei transformation'
1071	CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
1072	ELSEIF ( use_ug_for_galilei_tr .AND. &
1073	vg_vertical_gradient_level(1) /= 0.0 ) THEN
1074	message_string = 'baroclinicity (vg) not allowed simultaneously' // &
1075	' with galilei transformation'
1076	CALL message( 'check_parameters', 'PA0047', 1, 2, 0, 6, 0 )
1077	ELSE
1078	message_string = 'variable translation speed used for galilei-' // &
1079	'transformation, which may cause & instabilities in stably ' // &
1080	'stratified regions'
1081	CALL message( 'check_parameters', 'PA0048', 0, 1, 0, 6, 0 )
1082	ENDIF
1083	ENDIF
1084
1085	!
1086	!-- In case of using a prandtl-layer, calculated (or prescribed) surface
1087	!-- fluxes have to be used in the diffusion-terms
1088	IF ( prandtl_layer ) use_surface_fluxes = .TRUE.
1089
1090	!
1091	!-- Check boundary conditions and set internal variables:
1092	!-- Lateral boundary conditions
1093	IF ( bc_lr /= 'cyclic' .AND. bc_lr /= 'dirichlet/radiation' .AND. &
1094	bc_lr /= 'radiation/dirichlet' ) THEN
1095	message_string = 'unknown boundary condition: bc_lr = "' // &
1096	TRIM( bc_lr ) // '"'
1097	CALL message( 'check_parameters', 'PA0049', 1, 2, 0, 6, 0 )
1098	ENDIF
1099	IF ( bc_ns /= 'cyclic' .AND. bc_ns /= 'dirichlet/radiation' .AND. &
1100	bc_ns /= 'radiation/dirichlet' ) THEN
1101	message_string = 'unknown boundary condition: bc_ns = "' // &
1102	TRIM( bc_ns ) // '"'
1103	CALL message( 'check_parameters', 'PA0050', 1, 2, 0, 6, 0 )
1104	ENDIF
1105
1106	!
1107	!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
1108	!-- Willimas advection scheme. Several schemes and tools do not work with
1109	!-- non-cyclic boundary conditions.
1110	IF ( bc_lr /= 'cyclic' .OR. bc_ns /= 'cyclic' ) THEN
1111	IF ( psolver /= 'multigrid' ) THEN
1112	message_string = 'non-cyclic lateral boundaries do not allow ' // &
1113	'psolver = "' // TRIM( psolver ) // '"'
1114	CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
1115	ENDIF
1116	IF ( momentum_advec /= 'pw-scheme' ) THEN
1117	message_string = 'non-cyclic lateral boundaries do not allow ' // &
1118	'momentum_advec = "' // TRIM( momentum_advec ) // '"'
1119	CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
1120	ENDIF
1121	IF ( scalar_advec /= 'pw-scheme' ) THEN
1122	message_string = 'non-cyclic lateral boundaries do not allow ' // &
1123	'scalar_advec = "' // TRIM( scalar_advec ) // '"'
1124	CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
1125	ENDIF
1126	IF ( galilei_transformation ) THEN
1127	message_string = 'non-cyclic lateral boundaries do not allow ' // &
1128	'galilei_transformation = .T.'
1129	CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 )
1130	ENDIF
1131	ENDIF
1132
1133	!
1134	!-- Bottom boundary condition for the turbulent Kinetic energy
1135	IF ( bc_e_b == 'neumann' ) THEN
1136	ibc_e_b = 1
1137	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
1138	message_string = 'adjust_mixing_length = TRUE and bc_e_b = "neumann"'
1139	CALL message( 'check_parameters', 'PA0055', 0, 1, 0, 6, 0 )
1140	ENDIF
1141	ELSEIF ( bc_e_b == '(u)*2+neumann' ) THEN
1142	ibc_e_b = 2
1143	IF ( .NOT. adjust_mixing_length .AND. prandtl_layer ) THEN
1144	message_string = 'adjust_mixing_length = FALSE and bc_e_b = "' // &
1145	TRIM( bc_e_b ) // '"'
1146	CALL message( 'check_parameters', 'PA0056', 0, 1, 0, 6, 0 )
1147	ENDIF
1148	IF ( .NOT. prandtl_layer ) THEN
1149	bc_e_b = 'neumann'
1150	ibc_e_b = 1
1151	message_string = 'boundary condition bc_e_b changed to "' // &
1152	TRIM( bc_e_b ) // '"'
1153	CALL message( 'check_parameters', 'PA0057', 0, 1, 0, 6, 0 )
1154	ENDIF
1155	ELSE
1156	message_string = 'unknown boundary condition: bc_e_b = "' // &
1157	TRIM( bc_e_b ) // '"'
1158	CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 )
1159	ENDIF
1160
1161	!
1162	!-- Boundary conditions for perturbation pressure
1163	IF ( bc_p_b == 'dirichlet' ) THEN
1164	ibc_p_b = 0
1165	ELSEIF ( bc_p_b == 'neumann' ) THEN
1166	ibc_p_b = 1
1167	ELSEIF ( bc_p_b == 'neumann+inhomo' ) THEN
1168	ibc_p_b = 2
1169	ELSE
1170	message_string = 'unknown boundary condition: bc_p_b = "' // &
1171	TRIM( bc_p_b ) // '"'
1172	CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
1173	ENDIF
1174	IF ( ibc_p_b == 2 .AND. .NOT. prandtl_layer ) THEN
1175	message_string = 'boundary condition: bc_p_b = "' // TRIM( bc_p_b ) // &
1176	'" not allowed with prandtl_layer = .FALSE.'
1177	CALL message( 'check_parameters', 'PA0060', 1, 2, 0, 6, 0 )
1178	ENDIF
1179	IF ( bc_p_t == 'dirichlet' ) THEN
1180	ibc_p_t = 0
1181	ELSEIF ( bc_p_t == 'neumann' ) THEN
1182	ibc_p_t = 1
1183	ELSE
1184	message_string = 'unknown boundary condition: bc_p_t = "' // &
1185	TRIM( bc_p_t ) // '"'
1186	CALL message( 'check_parameters', 'PA0061', 1, 2, 0, 6, 0 )
1187	ENDIF
1188
1189	!
1190	!-- Boundary conditions for potential temperature
1191	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
1192	ibc_pt_b = 2
1193	ELSE
1194	IF ( bc_pt_b == 'dirichlet' ) THEN
1195	ibc_pt_b = 0
1196	ELSEIF ( bc_pt_b == 'neumann' ) THEN
1197	ibc_pt_b = 1
1198	ELSE
1199	message_string = 'unknown boundary condition: bc_pt_b = "' // &
1200	TRIM( bc_pt_b ) // '"'
1201	CALL message( 'check_parameters', 'PA0062', 1, 2, 0, 6, 0 )
1202	ENDIF
1203	ENDIF
1204
1205	IF ( bc_pt_t == 'dirichlet' ) THEN
1206	ibc_pt_t = 0
1207	ELSEIF ( bc_pt_t == 'neumann' ) THEN
1208	ibc_pt_t = 1
1209	ELSEIF ( bc_pt_t == 'initial_gradient' ) THEN
1210	ibc_pt_t = 2
1211	ELSE
1212	message_string = 'unknown boundary condition: bc_pt_t = "' // &
1213	TRIM( bc_pt_t ) // '"'
1214	CALL message( 'check_parameters', 'PA0063', 1, 2, 0, 6, 0 )
1215	ENDIF
1216
1217	IF ( surface_heatflux == 9999999.9 ) constant_heatflux = .FALSE.
1218	IF ( top_heatflux == 9999999.9 ) constant_top_heatflux = .FALSE.
1219	IF ( top_momentumflux_u /= 9999999.9 .AND. &
1220	top_momentumflux_v /= 9999999.9 ) THEN
1221	constant_top_momentumflux = .TRUE.
1222	ELSEIF ( .NOT. ( top_momentumflux_u == 9999999.9 .AND. &
1223	top_momentumflux_v == 9999999.9 ) ) THEN
1224	message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' // &
1225	'must be set'
1226	CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 )
1227	ENDIF
1228
1229	!
1230	!-- A given surface temperature implies Dirichlet boundary condition for
1231	!-- temperature. In this case specification of a constant heat flux is
1232	!-- forbidden.
1233	IF ( ibc_pt_b == 0 .AND. constant_heatflux .AND. &
1234	surface_heatflux /= 0.0 ) THEN
1235	message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //&
1236	'& is not allowed with constant_heatflux = .TRUE.'
1237	CALL message( 'check_parameters', 'PA0065', 1, 2, 0, 6, 0 )
1238	ENDIF
1239	IF ( constant_heatflux .AND. pt_surface_initial_change /= 0.0 ) THEN
1240	WRITE ( message_string, * ) 'constant_heatflux = .TRUE. is not allo', &
1241	'wed with pt_surface_initial_change (/=0) = ', &
1242	pt_surface_initial_change
1243	CALL message( 'check_parameters', 'PA0066', 1, 2, 0, 6, 0 )
1244	ENDIF
1245
1246	!
1247	!-- A given temperature at the top implies Dirichlet boundary condition for
1248	!-- temperature. In this case specification of a constant heat flux is
1249	!-- forbidden.
1250	IF ( ibc_pt_t == 0 .AND. constant_top_heatflux .AND. &
1251	top_heatflux /= 0.0 ) THEN
1252	message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //&
1253	'" is not allowed with constant_top_heatflux = .TRUE.'
1254	CALL message( 'check_parameters', 'PA0067', 1, 2, 0, 6, 0 )
1255	ENDIF
1256
1257	!
1258	!-- Boundary conditions for salinity
1259	IF ( ocean ) THEN
1260	IF ( bc_sa_t == 'dirichlet' ) THEN
1261	ibc_sa_t = 0
1262	ELSEIF ( bc_sa_t == 'neumann' ) THEN
1263	ibc_sa_t = 1
1264	ELSE
1265	message_string = 'unknown boundary condition: bc_sa_t = "' // &
1266	TRIM( bc_sa_t ) // '"'
1267	CALL message( 'check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
1268	ENDIF
1269
1270	IF ( top_salinityflux == 9999999.9 ) constant_top_salinityflux = .FALSE.
1271	IF ( ibc_sa_t == 1 .AND. top_salinityflux == 9999999.9 ) THEN
1272	message_string = 'boundary condition: bc_sa_t = "' // &
1273	TRIM( bc_sa_t ) // '" requires to set ' // &
1274	'top_salinityflux'
1275	CALL message( 'check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
1276	ENDIF
1277
1278	!
1279	!-- A fixed salinity at the top implies Dirichlet boundary condition for
1280	!-- salinity. In this case specification of a constant salinity flux is
1281	!-- forbidden.
1282	IF ( ibc_sa_t == 0 .AND. constant_top_salinityflux .AND. &
1283	top_salinityflux /= 0.0 ) THEN
1284	message_string = 'boundary condition: bc_sa_t = "' // &
1285	TRIM( bc_sa_t ) // '" is not allowed with ' // &
1286	'constant_top_salinityflux = .TRUE.'
1287	CALL message( 'check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
1288	ENDIF
1289
1290	ENDIF
1291
1292	!
1293	!-- In case of humidity or passive scalar, set boundary conditions for total
1294	!-- water content / scalar
1295	IF ( humidity .OR. passive_scalar ) THEN
1296	IF ( humidity ) THEN
1297	sq = 'q'
1298	ELSE
1299	sq = 's'
1300	ENDIF
1301	IF ( bc_q_b == 'dirichlet' ) THEN
1302	ibc_q_b = 0
1303	ELSEIF ( bc_q_b == 'neumann' ) THEN
1304	ibc_q_b = 1
1305	ELSE
1306	message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &
1307	'_b ="' // TRIM( bc_q_b ) // '"'
1308	CALL message( 'check_parameters', 'PA0071', 1, 2, 0, 6, 0 )
1309	ENDIF
1310	IF ( bc_q_t == 'dirichlet' ) THEN
1311	ibc_q_t = 0
1312	ELSEIF ( bc_q_t == 'neumann' ) THEN
1313	ibc_q_t = 1
1314	ELSE
1315	message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &
1316	'_t ="' // TRIM( bc_q_t ) // '"'
1317	CALL message( 'check_parameters', 'PA0072', 1, 2, 0, 6, 0 )
1318	ENDIF
1319
1320	IF ( surface_waterflux == 0.0 ) constant_waterflux = .FALSE.
1321
1322	!
1323	!-- A given surface humidity implies Dirichlet boundary condition for
1324	!-- humidity. In this case specification of a constant water flux is
1325	!-- forbidden.
1326	IF ( ibc_q_b == 0 .AND. constant_waterflux ) THEN
1327	message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b ' // &
1328	'= "' // TRIM( bc_q_b ) // '" is not allowed wi' // &
1329	'th prescribed surface flux'
1330	CALL message( 'check_parameters', 'PA0073', 1, 2, 0, 6, 0 )
1331	ENDIF
1332	IF ( constant_waterflux .AND. q_surface_initial_change /= 0.0 ) THEN
1333	WRITE( message_string, * ) 'a prescribed surface flux is not allo', &
1334	'wed with ', sq, '_surface_initial_change (/=0) = ', &
1335	q_surface_initial_change
1336	CALL message( 'check_parameters', 'PA0074', 1, 2, 0, 6, 0 )
1337	ENDIF
1338
1339	ENDIF
1340
1341	!
1342	!-- Boundary conditions for horizontal components of wind speed
1343	IF ( bc_uv_b == 'dirichlet' ) THEN
1344	ibc_uv_b = 0
1345	ELSEIF ( bc_uv_b == 'neumann' ) THEN
1346	ibc_uv_b = 1
1347	IF ( prandtl_layer ) THEN
1348	message_string = 'boundary condition: bc_uv_b = "' // &
1349	TRIM( bc_uv_b ) // '" is not allowed with prandtl_layer = .TRUE.'
1350	CALL message( 'check_parameters', 'PA0075', 1, 2, 0, 6, 0 )
1351	ENDIF
1352	ELSE
1353	message_string = 'unknown boundary condition: bc_uv_b = "' // &
1354	TRIM( bc_uv_b ) // '"'
1355	CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
1356	ENDIF
1357
1358	IF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
1359	bc_uv_t = 'neumann'
1360	ibc_uv_t = 1
1361	ELSE
1362	IF ( bc_uv_t == 'dirichlet' .OR. bc_uv_t == 'dirichlet_0' ) THEN
1363	ibc_uv_t = 0
1364	ELSEIF ( bc_uv_t == 'neumann' ) THEN
1365	ibc_uv_t = 1
1366	ELSE
1367	message_string = 'unknown boundary condition: bc_uv_t = "' // &
1368	TRIM( bc_uv_t ) // '"'
1369	CALL message( 'check_parameters', 'PA0077', 1, 2, 0, 6, 0 )
1370	ENDIF
1371	ENDIF
1372
1373	!
1374	!-- Compute and check, respectively, the Rayleigh Damping parameter
1375	IF ( rayleigh_damping_factor == -1.0 ) THEN
1376	IF ( momentum_advec == 'ups-scheme' ) THEN
1377	rayleigh_damping_factor = 0.01
1378	ELSE
1379	rayleigh_damping_factor = 0.0
1380	ENDIF
1381	ELSE
1382	IF ( rayleigh_damping_factor < 0.0 .OR. rayleigh_damping_factor > 1.0 ) &
1383	THEN
1384	WRITE( message_string, * ) 'rayleigh_damping_factor = ', &
1385	rayleigh_damping_factor, ' out of range [0.0,1.0]'
1386	CALL message( 'check_parameters', 'PA0078', 1, 2, 0, 6, 0 )
1387	ENDIF
1388	ENDIF
1389
1390	IF ( rayleigh_damping_height == -1.0 ) THEN
1391	IF ( .NOT. ocean ) THEN
1392	rayleigh_damping_height = 0.66666666666 * zu(nzt)
1393	ELSE
1394	rayleigh_damping_height = 0.66666666666 * zu(nzb)
1395	ENDIF
1396	ELSE
1397	IF ( .NOT. ocean ) THEN
1398	IF ( rayleigh_damping_height < 0.0 .OR. &
1399	rayleigh_damping_height > zu(nzt) ) THEN
1400	WRITE( message_string, * ) 'rayleigh_damping_height = ', &
1401	rayleigh_damping_height, ' out of range [0.0,', zu(nzt), ']'
1402	CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
1403	ENDIF
1404	ELSE
1405	IF ( rayleigh_damping_height > 0.0 .OR. &
1406	rayleigh_damping_height < zu(nzb) ) THEN
1407	WRITE( message_string, * ) 'rayleigh_damping_height = ', &
1408	rayleigh_damping_height, ' out of range [0.0,', zu(nzb), ']'
1409	CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
1410	ENDIF
1411	ENDIF
1412	ENDIF
1413
1414	!
1415	!-- Check limiters for Upstream-Spline scheme
1416	IF ( overshoot_limit_u < 0.0 .OR. overshoot_limit_v < 0.0 .OR. &
1417	overshoot_limit_w < 0.0 .OR. overshoot_limit_pt < 0.0 .OR. &
1418	overshoot_limit_e < 0.0 ) THEN
1419	message_string = 'overshoot_limit_... < 0.0 is not allowed'
1420	CALL message( 'check_parameters', 'PA0080', 1, 2, 0, 6, 0 )
1421	ENDIF
1422	IF ( ups_limit_u < 0.0 .OR. ups_limit_v < 0.0 .OR. ups_limit_w < 0.0 .OR. &
1423	ups_limit_pt < 0.0 .OR. ups_limit_e < 0.0 ) THEN
1424	message_string = 'ups_limit_... < 0.0 is not allowed'
1425	CALL message( 'check_parameters', 'PA0081', 1, 2, 0, 6, 0 )
1426	ENDIF
1427
1428	!
1429	!-- Check number of chosen statistic regions. More than 10 regions are not
1430	!-- allowed, because so far no more than 10 corresponding output files can
1431	!-- be opened (cf. check_open)
1432	IF ( statistic_regions > 9 .OR. statistic_regions < 0 ) THEN
1433	WRITE ( message_string, * ) 'number of statistic_regions = ', &
1434	statistic_regions+1, ' but only 10 regions are allowed'
1435	CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 )
1436	ENDIF
1437	IF ( normalizing_region > statistic_regions .OR. &
1438	normalizing_region < 0) THEN
1439	WRITE ( message_string, * ) 'normalizing_region = ', &
1440	normalizing_region, ' must be >= 0 and <= ',statistic_regions, &
1441	' (value of statistic_regions)'
1442	CALL message( 'check_parameters', 'PA0083', 1, 2, 0, 6, 0 )
1443	ENDIF
1444
1445	!
1446	!-- Check the interval for sorting particles.
1447	!-- Using particles as cloud droplets requires sorting after each timestep.
1448	IF ( dt_sort_particles /= 0.0 .AND. cloud_droplets ) THEN
1449	dt_sort_particles = 0.0
1450	message_string = 'dt_sort_particles is reset to 0.0 because of cloud' //&
1451	'_droplets = .TRUE.'
1452	CALL message( 'check_parameters', 'PA0084', 0, 1, 0, 6, 0 )
1453	ENDIF
1454
1455	!
1456	!-- Set the default intervals for data output, if necessary
1457	!-- NOTE: dt_dosp has already been set in package_parin
1458	IF ( dt_data_output /= 9999999.9 ) THEN
1459	IF ( dt_dopr == 9999999.9 ) dt_dopr = dt_data_output
1460	IF ( dt_dopts == 9999999.9 ) dt_dopts = dt_data_output
1461	IF ( dt_do2d_xy == 9999999.9 ) dt_do2d_xy = dt_data_output
1462	IF ( dt_do2d_xz == 9999999.9 ) dt_do2d_xz = dt_data_output
1463	IF ( dt_do2d_yz == 9999999.9 ) dt_do2d_yz = dt_data_output
1464	IF ( dt_do3d == 9999999.9 ) dt_do3d = dt_data_output
1465	IF ( dt_data_output_av == 9999999.9 ) dt_data_output_av = dt_data_output
1466	ENDIF
1467
1468	!
1469	!-- Set the default skip time intervals for data output, if necessary
1470	IF ( skip_time_dopr == 9999999.9 ) &
1471	skip_time_dopr = skip_time_data_output
1472	IF ( skip_time_dosp == 9999999.9 ) &
1473	skip_time_dosp = skip_time_data_output
1474	IF ( skip_time_do2d_xy == 9999999.9 ) &
1475	skip_time_do2d_xy = skip_time_data_output
1476	IF ( skip_time_do2d_xz == 9999999.9 ) &
1477	skip_time_do2d_xz = skip_time_data_output
1478	IF ( skip_time_do2d_yz == 9999999.9 ) &
1479	skip_time_do2d_yz = skip_time_data_output
1480	IF ( skip_time_do3d == 9999999.9 ) &
1481	skip_time_do3d = skip_time_data_output
1482	IF ( skip_time_data_output_av == 9999999.9 ) &
1483	skip_time_data_output_av = skip_time_data_output
1484
1485	!
1486	!-- Check the average intervals (first for 3d-data, then for profiles and
1487	!-- spectra)
1488	IF ( averaging_interval > dt_data_output_av ) THEN
1489	WRITE( message_string, * ) 'averaging_interval = ', &
1490	averaging_interval, ' must be <= dt_data_output = ', dt_data_output
1491	CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 )
1492	ENDIF
1493
1494	IF ( averaging_interval_pr == 9999999.9 ) THEN
1495	averaging_interval_pr = averaging_interval
1496	ENDIF
1497
1498	IF ( averaging_interval_pr > dt_dopr ) THEN
1499	WRITE( message_string, * ) 'averaging_interval_pr = ', &
1500	averaging_interval_pr, ' must be <= dt_dopr = ', dt_dopr
1501	CALL message( 'check_parameters', 'PA0086', 1, 2, 0, 6, 0 )
1502	ENDIF
1503
1504	IF ( averaging_interval_sp == 9999999.9 ) THEN
1505	averaging_interval_sp = averaging_interval
1506	ENDIF
1507
1508	IF ( averaging_interval_sp > dt_dosp ) THEN
1509	WRITE( message_string, * ) 'averaging_interval_sp = ', &
1510	averaging_interval_sp, ' must be <= dt_dosp = ', dt_dosp
1511	CALL message( 'check_parameters', 'PA0087', 1, 2, 0, 6, 0 )
1512	ENDIF
1513
1514	!
1515	!-- Set the default interval for profiles entering the temporal average
1516	IF ( dt_averaging_input_pr == 9999999.9 ) THEN
1517	dt_averaging_input_pr = dt_averaging_input
1518	ENDIF
1519
1520	!
1521	!-- Set the default interval for the output of timeseries to a reasonable
1522	!-- value (tries to minimize the number of calls of flow_statistics)
1523	IF ( dt_dots == 9999999.9 ) THEN
1524	IF ( averaging_interval_pr == 0.0 ) THEN
1525	dt_dots = MIN( dt_run_control, dt_dopr )
1526	ELSE
1527	dt_dots = MIN( dt_run_control, dt_averaging_input_pr )
1528	ENDIF
1529	ENDIF
1530
1531	!
1532	!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
1533	IF ( dt_averaging_input > averaging_interval ) THEN
1534	WRITE( message_string, * ) 'dt_averaging_input = ', &
1535	dt_averaging_input, ' must be <= averaging_interval = ', &
1536	averaging_interval
1537	CALL message( 'check_parameters', 'PA0088', 1, 2, 0, 6, 0 )
1538	ENDIF
1539
1540	IF ( dt_averaging_input_pr > averaging_interval_pr ) THEN
1541	WRITE( message_string, * ) 'dt_averaging_input_pr = ', &
1542	dt_averaging_input_pr, ' must be <= averaging_interval_pr = ', &
1543	averaging_interval_pr
1544	CALL message( 'check_parameters', 'PA0089', 1, 2, 0, 6, 0 )
1545	ENDIF
1546
1547	!
1548	!-- Set the default value for the integration interval of precipitation amount
1549	IF ( precipitation ) THEN
1550	IF ( precipitation_amount_interval == 9999999.9 ) THEN
1551	precipitation_amount_interval = dt_do2d_xy
1552	ELSE
1553	IF ( precipitation_amount_interval > dt_do2d_xy ) THEN
1554	WRITE( message_string, * ) 'precipitation_amount_interval = ', &
1555	precipitation_amount_interval, ' must not be larger than ', &
1556	'dt_do2d_xy = ', dt_do2d_xy
1557	CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
1558	ENDIF
1559	ENDIF
1560	ENDIF
1561
1562	!
1563	!-- Determine the number of output profiles and check whether they are
1564	!-- permissible
1565	DO WHILE ( data_output_pr(dopr_n+1) /= ' ' )
1566
1567	dopr_n = dopr_n + 1
1568	i = dopr_n
1569
1570	!
1571	!-- Determine internal profile number (for hom, homs)
1572	!-- and store height levels
1573	SELECT CASE ( TRIM( data_output_pr(i) ) )
1574
1575	CASE ( 'u', '#u' )
1576	dopr_index(i) = 1
1577	dopr_unit(i) = 'm/s'
1578	hom(:,2,1,:) = SPREAD( zu, 2, statistic_regions+1 )
1579	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1580	dopr_initial_index(i) = 5
1581	hom(:,2,5,:) = SPREAD( zu, 2, statistic_regions+1 )
1582	data_output_pr(i) = data_output_pr(i)(2:)
1583	ENDIF
1584
1585	CASE ( 'v', '#v' )
1586	dopr_index(i) = 2
1587	dopr_unit(i) = 'm/s'
1588	hom(:,2,2,:) = SPREAD( zu, 2, statistic_regions+1 )
1589	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1590	dopr_initial_index(i) = 6
1591	hom(:,2,6,:) = SPREAD( zu, 2, statistic_regions+1 )
1592	data_output_pr(i) = data_output_pr(i)(2:)
1593	ENDIF
1594
1595	CASE ( 'w' )
1596	dopr_index(i) = 3
1597	dopr_unit(i) = 'm/s'
1598	hom(:,2,3,:) = SPREAD( zw, 2, statistic_regions+1 )
1599
1600	CASE ( 'pt', '#pt' )
1601	IF ( .NOT. cloud_physics ) THEN
1602	dopr_index(i) = 4
1603	dopr_unit(i) = 'K'
1604	hom(:,2,4,:) = SPREAD( zu, 2, statistic_regions+1 )
1605	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1606	dopr_initial_index(i) = 7
1607	hom(:,2,7,:) = SPREAD( zu, 2, statistic_regions+1 )
1608	hom(nzb,2,7,:) = 0.0 ! because zu(nzb) is negative
1609	data_output_pr(i) = data_output_pr(i)(2:)
1610	ENDIF
1611	ELSE
1612	dopr_index(i) = 43
1613	dopr_unit(i) = 'K'
1614	hom(:,2,43,:) = SPREAD( zu, 2, statistic_regions+1 )
1615	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1616	dopr_initial_index(i) = 28
1617	hom(:,2,28,:) = SPREAD( zu, 2, statistic_regions+1 )
1618	hom(nzb,2,28,:) = 0.0 ! because zu(nzb) is negative
1619	data_output_pr(i) = data_output_pr(i)(2:)
1620	ENDIF
1621	ENDIF
1622
1623	CASE ( 'e' )
1624	dopr_index(i) = 8
1625	dopr_unit(i) = 'm2/s2'
1626	hom(:,2,8,:) = SPREAD( zu, 2, statistic_regions+1 )
1627	hom(nzb,2,8,:) = 0.0
1628
1629	CASE ( 'km', '#km' )
1630	dopr_index(i) = 9
1631	dopr_unit(i) = 'm2/s'
1632	hom(:,2,9,:) = SPREAD( zu, 2, statistic_regions+1 )
1633	hom(nzb,2,9,:) = 0.0
1634	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1635	dopr_initial_index(i) = 23
1636	hom(:,2,23,:) = hom(:,2,9,:)
1637	data_output_pr(i) = data_output_pr(i)(2:)
1638	ENDIF
1639
1640	CASE ( 'kh', '#kh' )
1641	dopr_index(i) = 10
1642	dopr_unit(i) = 'm2/s'
1643	hom(:,2,10,:) = SPREAD( zu, 2, statistic_regions+1 )
1644	hom(nzb,2,10,:) = 0.0
1645	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1646	dopr_initial_index(i) = 24
1647	hom(:,2,24,:) = hom(:,2,10,:)
1648	data_output_pr(i) = data_output_pr(i)(2:)
1649	ENDIF
1650
1651	CASE ( 'l', '#l' )
1652	dopr_index(i) = 11
1653	dopr_unit(i) = 'm'
1654	hom(:,2,11,:) = SPREAD( zu, 2, statistic_regions+1 )
1655	hom(nzb,2,11,:) = 0.0
1656	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1657	dopr_initial_index(i) = 25
1658	hom(:,2,25,:) = hom(:,2,11,:)
1659	data_output_pr(i) = data_output_pr(i)(2:)
1660	ENDIF
1661
1662	CASE ( 'w"u"' )
1663	dopr_index(i) = 12
1664	dopr_unit(i) = 'm2/s2'
1665	hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
1666	IF ( prandtl_layer ) hom(nzb,2,12,:) = zu(1)
1667
1668	CASE ( 'wu' )
1669	dopr_index(i) = 13
1670	dopr_unit(i) = 'm2/s2'
1671	hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
1672
1673	CASE ( 'w"v"' )
1674	dopr_index(i) = 14
1675	dopr_unit(i) = 'm2/s2'
1676	hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
1677	IF ( prandtl_layer ) hom(nzb,2,14,:) = zu(1)
1678
1679	CASE ( 'wv' )
1680	dopr_index(i) = 15
1681	dopr_unit(i) = 'm2/s2'
1682	hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
1683
1684	CASE ( 'w"pt"' )
1685	dopr_index(i) = 16
1686	dopr_unit(i) = 'K m/s'
1687	hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
1688
1689	CASE ( 'wpt' )
1690	dopr_index(i) = 17
1691	dopr_unit(i) = 'K m/s'
1692	hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
1693
1694	CASE ( 'wpt' )
1695	dopr_index(i) = 18
1696	dopr_unit(i) = 'K m/s'
1697	hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
1698
1699	CASE ( 'wu' )
1700	dopr_index(i) = 19
1701	dopr_unit(i) = 'm2/s2'
1702	hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
1703	IF ( prandtl_layer ) hom(nzb,2,19,:) = zu(1)
1704
1705	CASE ( 'wv' )
1706	dopr_index(i) = 20
1707	dopr_unit(i) = 'm2/s2'
1708	hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
1709	IF ( prandtl_layer ) hom(nzb,2,20,:) = zu(1)
1710
1711	CASE ( 'wptBC' )
1712	dopr_index(i) = 21
1713	dopr_unit(i) = 'K m/s'
1714	hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
1715
1716	CASE ( 'wptBC' )
1717	dopr_index(i) = 22
1718	dopr_unit(i) = 'K m/s'
1719	hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
1720
1721	CASE ( 'sa', '#sa' )
1722	IF ( .NOT. ocean ) THEN
1723	message_string = 'data_output_pr = ' // &
1724	TRIM( data_output_pr(i) ) // ' is not imp' // &
1725	'lemented for ocean = .FALSE.'
1726	CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
1727	ELSE
1728	dopr_index(i) = 23
1729	dopr_unit(i) = 'psu'
1730	hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 )
1731	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1732	dopr_initial_index(i) = 26
1733	hom(:,2,26,:) = SPREAD( zu, 2, statistic_regions+1 )
1734	hom(nzb,2,26,:) = 0.0 ! weil zu(nzb) negativ ist
1735	data_output_pr(i) = data_output_pr(i)(2:)
1736	ENDIF
1737	ENDIF
1738
1739	CASE ( 'u*2' )
1740	dopr_index(i) = 30
1741	dopr_unit(i) = 'm2/s2'
1742	hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
1743
1744	CASE ( 'v*2' )
1745	dopr_index(i) = 31
1746	dopr_unit(i) = 'm2/s2'
1747	hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
1748
1749	CASE ( 'w*2' )
1750	dopr_index(i) = 32
1751	dopr_unit(i) = 'm2/s2'
1752	hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
1753
1754	CASE ( 'pt*2' )
1755	dopr_index(i) = 33
1756	dopr_unit(i) = 'K2'
1757	hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
1758
1759	CASE ( 'e*' )
1760	dopr_index(i) = 34
1761	dopr_unit(i) = 'm2/s2'
1762	hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
1763
1764	CASE ( 'w2pt' )
1765	dopr_index(i) = 35
1766	dopr_unit(i) = 'K m2/s2'
1767	hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
1768
1769	CASE ( 'wpt2' )
1770	dopr_index(i) = 36
1771	dopr_unit(i) = 'K2 m/s'
1772	hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
1773
1774	CASE ( 'we' )
1775	dopr_index(i) = 37
1776	dopr_unit(i) = 'm3/s3'
1777	hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
1778
1779	CASE ( 'w*3' )
1780	dopr_index(i) = 38
1781	dopr_unit(i) = 'm3/s3'
1782	hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
1783
1784	CASE ( 'Sw' )
1785	dopr_index(i) = 39
1786	dopr_unit(i) = 'none'
1787	hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
1788
1789	CASE ( 'p' )
1790	dopr_index(i) = 40
1791	dopr_unit(i) = 'Pa'
1792	hom(:,2,40,:) = SPREAD( zu, 2, statistic_regions+1 )
1793
1794	CASE ( 'q', '#q' )
1795	IF ( .NOT. humidity ) THEN
1796	message_string = 'data_output_pr = ' // &
1797	TRIM( data_output_pr(i) ) // ' is not imp' // &
1798	'lemented for humidity = .FALSE.'
1799	CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
1800	ELSE
1801	dopr_index(i) = 41
1802	dopr_unit(i) = 'kg/kg'
1803	hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
1804	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1805	dopr_initial_index(i) = 26
1806	hom(:,2,26,:) = SPREAD( zu, 2, statistic_regions+1 )
1807	hom(nzb,2,26,:) = 0.0 ! weil zu(nzb) negativ ist
1808	data_output_pr(i) = data_output_pr(i)(2:)
1809	ENDIF
1810	ENDIF
1811
1812	CASE ( 's', '#s' )
1813	IF ( .NOT. passive_scalar ) THEN
1814	message_string = 'data_output_pr = ' // &
1815	TRIM( data_output_pr(i) ) // ' is not imp' // &
1816	'lemented for passive_scalar = .FALSE.'
1817	CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
1818	ELSE
1819	dopr_index(i) = 41
1820	dopr_unit(i) = 'kg/m3'
1821	hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
1822	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1823	dopr_initial_index(i) = 26
1824	hom(:,2,26,:) = SPREAD( zu, 2, statistic_regions+1 )
1825	hom(nzb,2,26,:) = 0.0 ! weil zu(nzb) negativ ist
1826	data_output_pr(i) = data_output_pr(i)(2:)
1827	ENDIF
1828	ENDIF
1829
1830	CASE ( 'qv', '#qv' )
1831	IF ( .NOT. cloud_physics ) THEN
1832	dopr_index(i) = 41
1833	dopr_unit(i) = 'kg/kg'
1834	hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
1835	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1836	dopr_initial_index(i) = 26
1837	hom(:,2,26,:) = SPREAD( zu, 2, statistic_regions+1 )
1838	hom(nzb,2,26,:) = 0.0 ! weil zu(nzb) negativ ist
1839	data_output_pr(i) = data_output_pr(i)(2:)
1840	ENDIF
1841	ELSE
1842	dopr_index(i) = 42
1843	dopr_unit(i) = 'kg/kg'
1844	hom(:,2,42,:) = SPREAD( zu, 2, statistic_regions+1 )
1845	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1846	dopr_initial_index(i) = 27
1847	hom(:,2,27,:) = SPREAD( zu, 2, statistic_regions+1 )
1848	hom(nzb,2,27,:) = 0.0 ! weil zu(nzb) negativ ist
1849	data_output_pr(i) = data_output_pr(i)(2:)
1850	ENDIF
1851	ENDIF
1852
1853	CASE ( 'lpt', '#lpt' )
1854	IF ( .NOT. cloud_physics ) THEN
1855	message_string = 'data_output_pr = ' // &
1856	TRIM( data_output_pr(i) ) // ' is not imp' // &
1857	'lemented for cloud_physics = .FALSE.'
1858	CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
1859	ELSE
1860	dopr_index(i) = 4
1861	dopr_unit(i) = 'K'
1862	hom(:,2,4,:) = SPREAD( zu, 2, statistic_regions+1 )
1863	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1864	dopr_initial_index(i) = 7
1865	hom(:,2,7,:) = SPREAD( zu, 2, statistic_regions+1 )
1866	hom(nzb,2,7,:) = 0.0 ! weil zu(nzb) negativ ist
1867	data_output_pr(i) = data_output_pr(i)(2:)
1868	ENDIF
1869	ENDIF
1870
1871	CASE ( 'vpt', '#vpt' )
1872	dopr_index(i) = 44
1873	dopr_unit(i) = 'K'
1874	hom(:,2,44,:) = SPREAD( zu, 2, statistic_regions+1 )
1875	IF ( data_output_pr(i)(1:1) == '#' ) THEN
1876	dopr_initial_index(i) = 29
1877	hom(:,2,29,:) = SPREAD( zu, 2, statistic_regions+1 )
1878	hom(nzb,2,29,:) = 0.0 ! weil zu(nzb) negativ ist
1879	data_output_pr(i) = data_output_pr(i)(2:)
1880	ENDIF
1881
1882	CASE ( 'w"vpt"' )
1883	dopr_index(i) = 45
1884	dopr_unit(i) = 'K m/s'
1885	hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
1886
1887	CASE ( 'wvpt' )
1888	dopr_index(i) = 46
1889	dopr_unit(i) = 'K m/s'
1890	hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
1891
1892	CASE ( 'wvpt' )
1893	dopr_index(i) = 47
1894	dopr_unit(i) = 'K m/s'
1895	hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
1896
1897	CASE ( 'w"q"' )
1898	IF ( .NOT. humidity ) THEN
1899	message_string = 'data_output_pr = ' // &
1900	TRIM( data_output_pr(i) ) // ' is not imp' // &
1901	'lemented for humidity = .FALSE.'
1902	CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
1903	ELSE
1904	dopr_index(i) = 48
1905	dopr_unit(i) = 'kg/kg m/s'
1906	hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
1907	ENDIF
1908
1909	CASE ( 'wq' )
1910	IF ( .NOT. humidity ) THEN
1911	message_string = 'data_output_pr = ' // &
1912	TRIM( data_output_pr(i) ) // ' is not imp' // &
1913	'lemented for humidity = .FALSE.'
1914	CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
1915	ELSE
1916	dopr_index(i) = 49
1917	dopr_unit(i) = 'kg/kg m/s'
1918	hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
1919	ENDIF
1920
1921	CASE ( 'wq' )
1922	IF ( .NOT. humidity ) THEN
1923	message_string = 'data_output_pr = ' // &
1924	TRIM( data_output_pr(i) ) // ' is not imp' // &
1925	'lemented for humidity = .FALSE.'
1926	CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
1927	ELSE
1928	dopr_index(i) = 50
1929	dopr_unit(i) = 'kg/kg m/s'
1930	hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
1931	ENDIF
1932
1933	CASE ( 'w"s"' )
1934	IF ( .NOT. passive_scalar ) THEN
1935	message_string = 'data_output_pr = ' // &
1936	TRIM( data_output_pr(i) ) // ' is not imp' // &
1937	'lemented for passive_scalar = .FALSE.'
1938	CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
1939	ELSE
1940	dopr_index(i) = 48
1941	dopr_unit(i) = 'kg/m3 m/s'
1942	hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
1943	ENDIF
1944
1945	CASE ( 'ws' )
1946	IF ( .NOT. passive_scalar ) THEN
1947	message_string = 'data_output_pr = ' // &
1948	TRIM( data_output_pr(i) ) // ' is not imp' // &
1949	'lemented for passive_scalar = .FALSE.'
1950	CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
1951	ELSE
1952	dopr_index(i) = 49
1953	dopr_unit(i) = 'kg/m3 m/s'
1954	hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
1955	ENDIF
1956
1957	CASE ( 'ws' )
1958	IF ( .NOT. passive_scalar ) THEN
1959	message_string = 'data_output_pr = ' // &
1960	TRIM( data_output_pr(i) ) // ' is not imp' // &
1961	'lemented for passive_scalar = .FALSE.'
1962	CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
1963	ELSE
1964	dopr_index(i) = 50
1965	dopr_unit(i) = 'kg/m3 m/s'
1966	hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
1967	ENDIF
1968
1969	CASE ( 'w"qv"' )
1970	IF ( humidity .AND. .NOT. cloud_physics ) &
1971	THEN
1972	dopr_index(i) = 48
1973	dopr_unit(i) = 'kg/kg m/s'
1974	hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
1975	ELSEIF( humidity .AND. cloud_physics ) THEN
1976	dopr_index(i) = 51
1977	dopr_unit(i) = 'kg/kg m/s'
1978	hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
1979	ELSE
1980	message_string = 'data_output_pr = ' // &
1981	TRIM( data_output_pr(i) ) // ' is not imp' // &
1982	'lemented for cloud_physics = .FALSE. an&' // &
1983	'd humidity = .FALSE.'
1984	CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
1985	ENDIF
1986
1987	CASE ( 'wqv' )
1988	IF ( humidity .AND. .NOT. cloud_physics ) &
1989	THEN
1990	dopr_index(i) = 49
1991	dopr_unit(i) = 'kg/kg m/s'
1992	hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
1993	ELSEIF( humidity .AND. cloud_physics ) THEN
1994	dopr_index(i) = 52
1995	dopr_unit(i) = 'kg/kg m/s'
1996	hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
1997	ELSE
1998	message_string = 'data_output_pr = ' // &
1999	TRIM( data_output_pr(i) ) // ' is not imp' // &
2000	'lemented for cloud_physics = .FALSE. an&' // &
2001	'd humidity = .FALSE.'
2002	CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2003	ENDIF
2004
2005	CASE ( 'wqv' )
2006	IF ( humidity .AND. .NOT. cloud_physics ) &
2007	THEN
2008	dopr_index(i) = 50
2009	dopr_unit(i) = 'kg/kg m/s'
2010	hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2011	ELSEIF( humidity .AND. cloud_physics ) THEN
2012	dopr_index(i) = 53
2013	dopr_unit(i) = 'kg/kg m/s'
2014	hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
2015	ELSE
2016	message_string = 'data_output_pr = ' // &
2017	TRIM( data_output_pr(i) ) // ' is not imp' // &
2018	'lemented for cloud_physics = .FALSE. an&' // &
2019	'd humidity = .FALSE.'
2020	CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2021	ENDIF
2022
2023	CASE ( 'ql' )
2024	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
2025	message_string = 'data_output_pr = ' // &
2026	TRIM( data_output_pr(i) ) // ' is not imp' // &
2027	'lemented for cloud_physics = .FALSE. or' // &
2028	'&cloud_droplets = .FALSE.'
2029	CALL message( 'check_parameters', 'PA0096', 1, 2, 0, 6, 0 )
2030	ELSE
2031	dopr_index(i) = 54
2032	dopr_unit(i) = 'kg/kg'
2033	hom(:,2,54,:) = SPREAD( zu, 2, statistic_regions+1 )
2034	ENDIF
2035
2036	CASE ( 'wuu*/dz' )
2037	dopr_index(i) = 55
2038	dopr_unit(i) = 'm2/s3'
2039	hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
2040
2041	CASE ( 'wp/dz' )
2042	dopr_index(i) = 56
2043	dopr_unit(i) = 'm2/s3'
2044	hom(:,2,56,:) = SPREAD( zw, 2, statistic_regions+1 )
2045
2046	CASE ( 'w"e/dz' )
2047	dopr_index(i) = 57
2048	dopr_unit(i) = 'm2/s3'
2049	hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
2050
2051	CASE ( 'u"pt"' )
2052	dopr_index(i) = 58
2053	dopr_unit(i) = 'K m/s'
2054	hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
2055
2056	CASE ( 'upt' )
2057	dopr_index(i) = 59
2058	dopr_unit(i) = 'K m/s'
2059	hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
2060
2061	CASE ( 'upt_t' )
2062	dopr_index(i) = 60
2063	dopr_unit(i) = 'K m/s'
2064	hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
2065
2066	CASE ( 'v"pt"' )
2067	dopr_index(i) = 61
2068	dopr_unit(i) = 'K m/s'
2069	hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
2070
2071	CASE ( 'vpt' )
2072	dopr_index(i) = 62
2073	dopr_unit(i) = 'K m/s'
2074	hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
2075
2076	CASE ( 'vpt_t' )
2077	dopr_index(i) = 63
2078	dopr_unit(i) = 'K m/s'
2079	hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
2080
2081	CASE ( 'rho' )
2082	dopr_index(i) = 64
2083	dopr_unit(i) = 'kg/m3'
2084	hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
2085
2086	CASE ( 'w"sa"' )
2087	IF ( .NOT. ocean ) THEN
2088	message_string = 'data_output_pr = ' // &
2089	TRIM( data_output_pr(i) ) // ' is not imp' // &
2090	'lemented for ocean = .FALSE.'
2091	CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2092	ELSE
2093	dopr_index(i) = 65
2094	dopr_unit(i) = 'psu m/s'
2095	hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 )
2096	ENDIF
2097
2098	CASE ( 'wsa' )
2099	IF ( .NOT. ocean ) THEN
2100	message_string = 'data_output_pr = ' // &
2101	TRIM( data_output_pr(i) ) // ' is not imp' // &
2102	'lemented for ocean = .FALSE.'
2103	CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2104	ELSE
2105	dopr_index(i) = 66
2106	dopr_unit(i) = 'psu m/s'
2107	hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 )
2108	ENDIF
2109
2110	CASE ( 'wsa' )
2111	IF ( .NOT. ocean ) THEN
2112	message_string = 'data_output_pr = ' // &
2113	TRIM( data_output_pr(i) ) // ' is not imp' // &
2114	'lemented for ocean = .FALSE.'
2115	CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2116	ELSE
2117	dopr_index(i) = 67
2118	dopr_unit(i) = 'psu m/s'
2119	hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
2120	ENDIF
2121
2122	CASE ( 'wp' )
2123	dopr_index(i) = 68
2124	dopr_unit(i) = 'm3/s3'
2125	hom(:,2,68,:) = SPREAD( zu, 2, statistic_regions+1 )
2126
2127	CASE ( 'w"e' )
2128	dopr_index(i) = 69
2129	dopr_unit(i) = 'm3/s3'
2130	hom(:,2,69,:) = SPREAD( zu, 2, statistic_regions+1 )
2131
2132	CASE ( 'q*2' )
2133	IF ( .NOT. humidity ) THEN
2134	message_string = 'data_output_pr = ' // &
2135	TRIM( data_output_pr(i) ) // ' is not imp' // &
2136	'lemented for humidity = .FALSE.'
2137	CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2138	ELSE
2139	dopr_index(i) = 70
2140	dopr_unit(i) = 'kg2/kg2'
2141	hom(:,2,70,:) = SPREAD( zu, 2, statistic_regions+1 )
2142	ENDIF
2143
2144	CASE DEFAULT
2145
2146	CALL user_check_data_output_pr( data_output_pr(i), i, unit )
2147
2148	IF ( unit == 'illegal' ) THEN
2149	IF ( data_output_pr_user(1) /= ' ' ) THEN
2150	message_string = 'illegal value for data_output_pr or ' // &
2151	'data_output_pr_user = "' // &
2152	TRIM( data_output_pr(i) ) // '"'
2153	CALL message( 'check_parameters', 'PA0097', 1, 2, 0, 6, 0 )
2154	ELSE
2155	message_string = 'illegal value for data_output_pr = "' // &
2156	TRIM( data_output_pr(i) ) // '"'
2157	CALL message( 'check_parameters', 'PA0098', 1, 2, 0, 6, 0 )
2158	ENDIF
2159	ENDIF
2160
2161	END SELECT
2162	!
2163	!-- Check to which of the predefined coordinate systems the profile belongs
2164	DO k = 1, crmax
2165	IF ( INDEX( cross_profiles(k), ' '//TRIM( data_output_pr(i) )//' ' ) &
2166	/=0 ) &
2167	THEN
2168	dopr_crossindex(i) = k
2169	EXIT
2170	ENDIF
2171	ENDDO
2172	!
2173	!-- Generate the text for the labels of the PROFIL output file. "-characters
2174	!-- must be substituted, otherwise PROFIL would interpret them as TeX
2175	!-- control characters
2176	dopr_label(i) = data_output_pr(i)
2177	position = INDEX( dopr_label(i) , '"' )
2178	DO WHILE ( position /= 0 )
2179	dopr_label(i)(position:position) = ''''
2180	position = INDEX( dopr_label(i) , '"' )
2181	ENDDO
2182
2183	ENDDO
2184
2185	!
2186	!-- y-value range of the coordinate system (PROFIL).
2187	!-- x-value range determined in plot_1d.
2188	IF ( .NOT. ocean ) THEN
2189	cross_uymin = 0.0
2190	IF ( z_max_do1d == -1.0 ) THEN
2191	cross_uymax = zu(nzt+1)
2192	ELSEIF ( z_max_do1d < zu(nzb+1) .OR. z_max_do1d > zu(nzt+1) ) THEN
2193	WRITE( message_string, * ) 'z_max_do1d = ', z_max_do1d, ' must ', &
2194	'be >= ', zu(nzb+1), ' or <= ', zu(nzt+1)
2195	CALL message( 'check_parameters', 'PA0099', 1, 2, 0, 6, 0 )
2196	ELSE
2197	cross_uymax = z_max_do1d
2198	ENDIF
2199	ENDIF
2200
2201	!
2202	!-- Check whether the chosen normalizing factor for the coordinate systems is
2203	!-- permissible
2204	DO i = 1, crmax
2205	SELECT CASE ( TRIM( cross_normalized_x(i) ) ) ! TRIM required on IBM
2206
2207	CASE ( '', 'wpt0', 'ws2', 'tsw2', 'ws3', 'ws2tsw', 'wstsw2' )
2208	j = 0
2209
2210	CASE DEFAULT
2211	message_string = 'unknown normalization method cross_normali' // &
2212	'zed_x = "' // TRIM( cross_normalized_x(i) ) // &
2213	'"'
2214	CALL message( 'check_parameters', 'PA0100', 1, 2, 0, 6, 0 )
2215
2216	END SELECT
2217	SELECT CASE ( TRIM( cross_normalized_y(i) ) ) ! TRIM required on IBM
2218
2219	CASE ( '', 'z_i' )
2220	j = 0
2221
2222	CASE DEFAULT
2223	message_string = 'unknown normalization method cross_normali' // &
2224	'zed_y = "' // TRIM( cross_normalized_y(i) ) // &
2225	'"'
2226	CALL message( 'check_parameters', 'PA0101', 1, 2, 0, 6, 0 )
2227
2228	END SELECT
2229	ENDDO
2230	!
2231	!-- Check normalized y-value range of the coordinate system (PROFIL)
2232	IF ( z_max_do1d_normalized /= -1.0 .AND. z_max_do1d_normalized <= 0.0 ) &
2233	THEN
2234	WRITE( message_string, * ) 'z_max_do1d_normalized = ', &
2235	z_max_do1d_normalized, ' must be >= 0.0'
2236	CALL message( 'check_parameters', 'PA0101', 1, 2, 0, 6, 0 )
2237	ENDIF
2238
2239
2240	!
2241	!-- Append user-defined data output variables to the standard data output
2242	IF ( data_output_user(1) /= ' ' ) THEN
2243	i = 1
2244	DO WHILE ( data_output(i) /= ' ' .AND. i <= 100 )
2245	i = i + 1
2246	ENDDO
2247	j = 1
2248	DO WHILE ( data_output_user(j) /= ' ' .AND. j <= 100 )
2249	IF ( i > 100 ) THEN
2250	message_string = 'number of output quantitities given by data' // &
2251	'_output and data_output_user exceeds the limit of 100'
2252	CALL message( 'check_parameters', 'PA0102', 1, 2, 0, 6, 0 )
2253	ENDIF
2254	data_output(i) = data_output_user(j)
2255	i = i + 1
2256	j = j + 1
2257	ENDDO
2258	ENDIF
2259
2260	!
2261	!-- Check and set steering parameters for 2d/3d data output and averaging
2262	i = 1
2263	DO WHILE ( data_output(i) /= ' ' .AND. i <= 100 )
2264	!
2265	!-- Check for data averaging
2266	ilen = LEN_TRIM( data_output(i) )
2267	j = 0 ! no data averaging
2268	IF ( ilen > 3 ) THEN
2269	IF ( data_output(i)(ilen-2:ilen) == '_av' ) THEN
2270	j = 1 ! data averaging
2271	data_output(i) = data_output(i)(1:ilen-3)
2272	ENDIF
2273	ENDIF
2274	!
2275	!-- Check for cross section or volume data
2276	ilen = LEN_TRIM( data_output(i) )
2277	k = 0 ! 3d data
2278	var = data_output(i)(1:ilen)
2279	IF ( ilen > 3 ) THEN
2280	IF ( data_output(i)(ilen-2:ilen) == '_xy' .OR. &
2281	data_output(i)(ilen-2:ilen) == '_xz' .OR. &
2282	data_output(i)(ilen-2:ilen) == '_yz' ) THEN
2283	k = 1 ! 2d data
2284	var = data_output(i)(1:ilen-3)
2285	ENDIF
2286	ENDIF
2287	!
2288	!-- Check for allowed value and set units
2289	SELECT CASE ( TRIM( var ) )
2290
2291	CASE ( 'e' )
2292	IF ( constant_diffusion ) THEN
2293	message_string = 'output of "' // TRIM( var ) // '" requi' // &
2294	'res constant_diffusion = .FALSE.'
2295	CALL message( 'check_parameters', 'PA0103', 1, 2, 0, 6, 0 )
2296	ENDIF
2297	unit = 'm2/s2'
2298
2299	CASE ( 'pc', 'pr' )
2300	IF ( .NOT. particle_advection ) THEN
2301	message_string = 'output of "' // TRIM( var ) // '" requir' // &
2302	'es a "particles_par"-NAMELIST in the parameter file (PARIN)'
2303	CALL message( 'check_parameters', 'PA0104', 1, 2, 0, 6, 0 )
2304	ENDIF
2305	IF ( TRIM( var ) == 'pc' ) unit = 'number'
2306	IF ( TRIM( var ) == 'pr' ) unit = 'm'
2307
2308	CASE ( 'q', 'vpt' )
2309	IF ( .NOT. humidity ) THEN
2310	message_string = 'output of "' // TRIM( var ) // '" requi' // &
2311	'res humidity = .TRUE.'
2312	CALL message( 'check_parameters', 'PA0105', 1, 2, 0, 6, 0 )
2313	ENDIF
2314	IF ( TRIM( var ) == 'q' ) unit = 'kg/kg'
2315	IF ( TRIM( var ) == 'vpt' ) unit = 'K'
2316
2317	CASE ( 'ql' )
2318	IF ( .NOT. ( cloud_physics .OR. cloud_droplets ) ) THEN
2319	message_string = 'output of "' // TRIM( var ) // '" requi' // &
2320	'res cloud_physics = .TRUE. or cloud_droplets = .TRUE.'
2321	CALL message( 'check_parameters', 'PA0106', 1, 2, 0, 6, 0 )
2322	ENDIF
2323	unit = 'kg/kg'
2324
2325	CASE ( 'ql_c', 'ql_v', 'ql_vp' )
2326	IF ( .NOT. cloud_droplets ) THEN
2327	message_string = 'output of "' // TRIM( var ) // '" requi' // &
2328	'res cloud_droplets = .TRUE.'
2329	CALL message( 'check_parameters', 'PA0107', 1, 2, 0, 6, 0 )
2330	ENDIF
2331	IF ( TRIM( var ) == 'ql_c' ) unit = 'kg/kg'
2332	IF ( TRIM( var ) == 'ql_v' ) unit = 'm3'
2333	IF ( TRIM( var ) == 'ql_vp' ) unit = 'none'
2334
2335	CASE ( 'qv' )
2336	IF ( .NOT. cloud_physics ) THEN
2337	message_string = 'output of "' // TRIM( var ) // '" requi' // &
2338	'res cloud_physics = .TRUE.'
2339	CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2340	ENDIF
2341	unit = 'kg/kg'
2342
2343	CASE ( 'rho' )
2344	IF ( .NOT. ocean ) THEN
2345	message_string = 'output of "' // TRIM( var ) // '" requi' // &
2346	'res ocean = .TRUE.'
2347	CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
2348	ENDIF
2349	unit = 'kg/m3'
2350
2351	CASE ( 's' )
2352	IF ( .NOT. passive_scalar ) THEN
2353	message_string = 'output of "' // TRIM( var ) // '" requi' // &
2354	'res passive_scalar = .TRUE.'
2355	CALL message( 'check_parameters', 'PA0110', 1, 2, 0, 6, 0 )
2356	ENDIF
2357	unit = 'conc'
2358
2359	CASE ( 'sa' )
2360	IF ( .NOT. ocean ) THEN
2361	message_string = 'output of "' // TRIM( var ) // '" requi' // &
2362	'res ocean = .TRUE.'
2363	CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
2364	ENDIF
2365	unit = 'psu'
2366
2367	CASE ( 'u', 't', 'lwp', 'pra', 'prr', 'qsws', 'shf', 'z0' )
2368	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
2369	message_string = 'illegal value for data_output: "' // &
2370	TRIM( var ) // '" & only 2d-horizontal ' // &
2371	'cross sections are allowed for this value'
2372	CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
2373	ENDIF
2374	IF ( TRIM( var ) == 'lwp*' .AND. .NOT. cloud_physics ) THEN
2375	message_string = 'output of "' // TRIM( var ) // '" requi' // &
2376	'res cloud_physics = .TRUE.'
2377	CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2378	ENDIF
2379	IF ( TRIM( var ) == 'pra*' .AND. .NOT. precipitation ) THEN
2380	message_string = 'output of "' // TRIM( var ) // '" requi' // &
2381	'res precipitation = .TRUE.'
2382	CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
2383	ENDIF
2384	IF ( TRIM( var ) == 'pra*' .AND. j == 1 ) THEN
2385	message_string = 'temporal averaging of precipitation ' // &
2386	'amount "' // TRIM( var ) // '" is not possible'
2387	CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
2388	ENDIF
2389	IF ( TRIM( var ) == 'prr*' .AND. .NOT. precipitation ) THEN
2390	message_string = 'output of "' // TRIM( var ) // '" requi' // &
2391	'res precipitation = .TRUE.'
2392	CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
2393	ENDIF
2394	IF ( TRIM( var ) == 'qsws*' .AND. .NOT. humidity ) THEN
2395	message_string = 'output of "' // TRIM( var ) // '" requi' // &
2396	'res humidity = .TRUE.'
2397	CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
2398	ENDIF
2399
2400	IF ( TRIM( var ) == 'lwp' ) unit = 'kg/kgm'
2401	IF ( TRIM( var ) == 'pra*' ) unit = 'mm'
2402	IF ( TRIM( var ) == 'prr*' ) unit = 'mm/s'
2403	IF ( TRIM( var ) == 'qsws*' ) unit = 'kgm/kgs'
2404	IF ( TRIM( var ) == 'shf' ) unit = 'Km/s'
2405	IF ( TRIM( var ) == 't*' ) unit = 'K'
2406	IF ( TRIM( var ) == 'u*' ) unit = 'm/s'
2407	IF ( TRIM( var ) == 'z0*' ) unit = 'm'
2408
2409
2410	CASE ( 'p', 'pt', 'u', 'v', 'w' )
2411	IF ( TRIM( var ) == 'p' ) unit = 'Pa'
2412	IF ( TRIM( var ) == 'pt' ) unit = 'K'
2413	IF ( TRIM( var ) == 'u' ) unit = 'm/s'
2414	IF ( TRIM( var ) == 'v' ) unit = 'm/s'
2415	IF ( TRIM( var ) == 'w' ) unit = 'm/s'
2416	CONTINUE
2417
2418	CASE DEFAULT
2419	CALL user_check_data_output( var, unit )
2420
2421	IF ( unit == 'illegal' ) THEN
2422	IF ( data_output_user(1) /= ' ' ) THEN
2423	message_string = 'illegal value for data_output or ' // &
2424	'data_output_user = "' // TRIM( data_output(i) ) // '"'
2425	CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
2426	ELSE
2427	message_string = 'illegal value for data_output =' // &
2428	TRIM( data_output(i) ) // '"'
2429	CALL message( 'check_parameters', 'PA0115', 1, 2, 0, 6, 0 )
2430	ENDIF
2431	ENDIF
2432
2433	END SELECT
2434	!
2435	!-- Set the internal steering parameters appropriately
2436	IF ( k == 0 ) THEN
2437	do3d_no(j) = do3d_no(j) + 1
2438	do3d(j,do3d_no(j)) = data_output(i)
2439	do3d_unit(j,do3d_no(j)) = unit
2440	ELSE
2441	do2d_no(j) = do2d_no(j) + 1
2442	do2d(j,do2d_no(j)) = data_output(i)
2443	do2d_unit(j,do2d_no(j)) = unit
2444	IF ( data_output(i)(ilen-2:ilen) == '_xy' ) THEN
2445	data_output_xy(j) = .TRUE.
2446	ENDIF
2447	IF ( data_output(i)(ilen-2:ilen) == '_xz' ) THEN
2448	data_output_xz(j) = .TRUE.
2449	ENDIF
2450	IF ( data_output(i)(ilen-2:ilen) == '_yz' ) THEN
2451	data_output_yz(j) = .TRUE.
2452	ENDIF
2453	ENDIF
2454
2455	IF ( j == 1 ) THEN
2456	!
2457	!-- Check, if variable is already subject to averaging
2458	found = .FALSE.
2459	DO k = 1, doav_n
2460	IF ( TRIM( doav(k) ) == TRIM( var ) ) found = .TRUE.
2461	ENDDO
2462
2463	IF ( .NOT. found ) THEN
2464	doav_n = doav_n + 1
2465	doav(doav_n) = var
2466	ENDIF
2467	ENDIF
2468
2469	i = i + 1
2470	ENDDO
2471
2472	!
2473	!-- Check sectional planes and store them in one shared array
2474	IF ( ANY( section_xy > nz + 1 ) ) THEN
2475	WRITE( message_string, * ) 'section_xy must be <= nz + 1 = ', nz + 1
2476	CALL message( 'check_parameters', 'PA0319', 1, 2, 0, 6, 0 )
2477	ENDIF
2478	IF ( ANY( section_xz > ny + 1 ) ) THEN
2479	WRITE( message_string, * ) 'section_xz must be <= ny + 1 = ', ny + 1
2480	CALL message( 'check_parameters', 'PA0320', 1, 2, 0, 6, 0 )
2481	ENDIF
2482	IF ( ANY( section_yz > nx + 1 ) ) THEN
2483	WRITE( message_string, * ) 'section_yz must be <= nx + 1 = ', nx + 1
2484	CALL message( 'check_parameters', 'PA0321', 1, 2, 0, 6, 0 )
2485	ENDIF
2486	section(:,1) = section_xy
2487	section(:,2) = section_xz
2488	section(:,3) = section_yz
2489
2490	!
2491	!-- Upper plot limit (grid point value) for 1D profiles
2492	IF ( z_max_do1d == -1.0 ) THEN
2493	nz_do1d = nzt+1
2494	ELSE
2495	DO k = nzb+1, nzt+1
2496	nz_do1d = k
2497	IF ( zw(k) > z_max_do1d ) EXIT
2498	ENDDO
2499	ENDIF
2500
2501	!
2502	!-- Upper plot limit for 2D vertical sections
2503	IF ( z_max_do2d == -1.0 ) z_max_do2d = zu(nzt)
2504	IF ( z_max_do2d < zu(nzb+1) .OR. z_max_do2d > zu(nzt) ) THEN
2505	WRITE( message_string, * ) 'z_max_do2d = ', z_max_do2d, &
2506	' must be >= ', zu(nzb+1), '(zu(nzb+1)) and <= ', zu(nzt), &
2507	' (zu(nzt))'
2508	CALL message( 'check_parameters', 'PA0116', 1, 2, 0, 6, 0 )
2509	ENDIF
2510
2511	!
2512	!-- Upper plot limit for 3D arrays
2513	IF ( nz_do3d == -9999 ) nz_do3d = nzt + 1
2514
2515	!
2516	!-- Determine and check accuracy for compressed 3D plot output
2517	IF ( do3d_compress ) THEN
2518	!
2519	!-- Compression only permissible on T3E machines
2520	IF ( host(1:3) /= 't3e' ) THEN
2521	message_string = 'do3d_compress = .TRUE. not allowed on host "' // &
2522	TRIM( host ) // '"'
2523	CALL message( 'check_parameters', 'PA0117', 1, 2, 0, 6, 0 )
2524	ENDIF
2525
2526	i = 1
2527	DO WHILE ( do3d_comp_prec(i) /= ' ' )
2528
2529	ilen = LEN_TRIM( do3d_comp_prec(i) )
2530	IF ( LLT( do3d_comp_prec(i)(ilen:ilen), '0' ) .OR. &
2531	LGT( do3d_comp_prec(i)(ilen:ilen), '9' ) ) THEN
2532	WRITE( message_string, * ) 'illegal precision: do3d_comp_prec', &
2533	'(', i, ') = "', TRIM(do3d_comp_prec(i)),'"'
2534	CALL message( 'check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
2535	ENDIF
2536
2537	prec = IACHAR( do3d_comp_prec(i)(ilen:ilen) ) - IACHAR( '0' )
2538	var = do3d_comp_prec(i)(1:ilen-1)
2539
2540	SELECT CASE ( var )
2541
2542	CASE ( 'u' )
2543	j = 1
2544	CASE ( 'v' )
2545	j = 2
2546	CASE ( 'w' )
2547	j = 3
2548	CASE ( 'p' )
2549	j = 4
2550	CASE ( 'pt' )
2551	j = 5
2552
2553	CASE DEFAULT
2554	WRITE( message_string, * ) 'unknown variable "', &
2555	TRIM( do3d_comp_prec(i) ), '" given for do3d_comp_prec(', &
2556	i, ')'
2557	CALL message( 'check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
2558
2559	END SELECT
2560
2561	plot_3d_precision(j)%precision = prec
2562	i = i + 1
2563
2564	ENDDO
2565	ENDIF
2566
2567	!
2568	!-- Check the data output format(s)
2569	IF ( data_output_format(1) == ' ' ) THEN
2570	!
2571	!-- Default value
2572	netcdf_output = .TRUE.
2573	ELSE
2574	i = 1
2575	DO WHILE ( data_output_format(i) /= ' ' )
2576
2577	SELECT CASE ( data_output_format(i) )
2578
2579	CASE ( 'netcdf' )
2580	netcdf_output = .TRUE.
2581	CASE ( 'iso2d' )
2582	iso2d_output = .TRUE.
2583	CASE ( 'profil' )
2584	profil_output = .TRUE.
2585	CASE ( 'avs' )
2586	avs_output = .TRUE.
2587
2588	CASE DEFAULT
2589	message_string = 'unknown value for data_output_format "' // &
2590	TRIM( data_output_format(i) ) // '"'
2591	CALL message( 'check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
2592
2593	END SELECT
2594
2595	i = i + 1
2596	IF ( i > 10 ) EXIT
2597
2598	ENDDO
2599
2600	ENDIF
2601
2602	!
2603	!-- Check netcdf precison
2604	ldum = .FALSE.
2605	CALL define_netcdf_header( 'ch', ldum, 0 )
2606
2607	!
2608	!-- Check, whether a constant diffusion coefficient shall be used
2609	IF ( km_constant /= -1.0 ) THEN
2610	IF ( km_constant < 0.0 ) THEN
2611	WRITE( message_string, * ) 'km_constant = ', km_constant, ' < 0.0'
2612	CALL message( 'check_parameters', 'PA0121', 1, 2, 0, 6, 0 )
2613	ELSE
2614	IF ( prandtl_number < 0.0 ) THEN
2615	WRITE( message_string, * ) 'prandtl_number = ', prandtl_number, &
2616	' < 0.0'
2617	CALL message( 'check_parameters', 'PA0122', 1, 2, 0, 6, 0 )
2618	ENDIF
2619	constant_diffusion = .TRUE.
2620
2621	IF ( prandtl_layer ) THEN
2622	message_string = 'prandtl_layer is not allowed with fixed ' // &
2623	'value of km'
2624	CALL message( 'check_parameters', 'PA0123', 1, 2, 0, 6, 0 )
2625	ENDIF
2626	ENDIF
2627	ENDIF
2628
2629	!
2630	!-- In case of non-cyclic lateral boundaries, set the default maximum value
2631	!-- for the horizontal diffusivity used within the outflow damping layer,
2632	!-- and check/set the width of the damping layer
2633	IF ( bc_lr /= 'cyclic' ) THEN
2634	IF ( km_damp_max == -1.0 ) THEN
2635	km_damp_max = 0.5 * dx
2636	ENDIF
2637	IF ( outflow_damping_width == -1.0 ) THEN
2638	outflow_damping_width = MIN( 20, nx/2 )
2639	ENDIF
2640	IF ( outflow_damping_width <= 0 .OR. outflow_damping_width > nx ) THEN
2641	message_string = 'outflow_damping width out of range'
2642	CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
2643	ENDIF
2644	ENDIF
2645
2646	IF ( bc_ns /= 'cyclic' ) THEN
2647	IF ( km_damp_max == -1.0 ) THEN
2648	km_damp_max = 0.5 * dy
2649	ENDIF
2650	IF ( outflow_damping_width == -1.0 ) THEN
2651	outflow_damping_width = MIN( 20, ny/2 )
2652	ENDIF
2653	IF ( outflow_damping_width <= 0 .OR. outflow_damping_width > ny ) THEN
2654	message_string = 'outflow_damping width out of range'
2655	CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
2656	ENDIF
2657	ENDIF
2658
2659	!
2660	!-- Check value range for rif
2661	IF ( rif_min >= rif_max ) THEN
2662	WRITE( message_string, * ) 'rif_min = ', rif_min, ' must be less ', &
2663	'than rif_max = ', rif_max
2664	CALL message( 'check_parameters', 'PA0125', 1, 2, 0, 6, 0 )
2665	ENDIF
2666
2667	!
2668	!-- Determine upper and lower hight level indices for random perturbations
2669	IF ( disturbance_level_b == -9999999.9 ) THEN
2670	IF ( ocean ) THEN
2671	disturbance_level_b = zu((nzt*2)/3)
2672	disturbance_level_ind_b = ( nzt * 2 ) / 3
2673	ELSE
2674	disturbance_level_b = zu(nzb+3)
2675	disturbance_level_ind_b = nzb + 3
2676	ENDIF
2677	ELSEIF ( disturbance_level_b < zu(3) ) THEN
2678	WRITE( message_string, * ) 'disturbance_level_b = ', &
2679	disturbance_level_b, ' must be >= ', zu(3), '(zu(3))'
2680	CALL message( 'check_parameters', 'PA0126', 1, 2, 0, 6, 0 )
2681	ELSEIF ( disturbance_level_b > zu(nzt-2) ) THEN
2682	WRITE( message_string, * ) 'disturbance_level_b = ', &
2683	disturbance_level_b, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
2684	CALL message( 'check_parameters', 'PA0127', 1, 2, 0, 6, 0 )
2685	ELSE
2686	DO k = 3, nzt-2
2687	IF ( disturbance_level_b <= zu(k) ) THEN
2688	disturbance_level_ind_b = k
2689	EXIT
2690	ENDIF
2691	ENDDO
2692	ENDIF
2693
2694	IF ( disturbance_level_t == -9999999.9 ) THEN
2695	IF ( ocean ) THEN
2696	disturbance_level_t = zu(nzt-3)
2697	disturbance_level_ind_t = nzt - 3
2698	ELSE
2699	disturbance_level_t = zu(nzt/3)
2700	disturbance_level_ind_t = nzt / 3
2701	ENDIF
2702	ELSEIF ( disturbance_level_t > zu(nzt-2) ) THEN
2703	WRITE( message_string, * ) 'disturbance_level_t = ', &
2704	disturbance_level_t, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
2705	CALL message( 'check_parameters', 'PA0128', 1, 2, 0, 6, 0 )
2706	ELSEIF ( disturbance_level_t < disturbance_level_b ) THEN
2707	WRITE( message_string, * ) 'disturbance_level_t = ', &
2708	disturbance_level_t, ' must be >= disturbance_level_b = ', &
2709	disturbance_level_b
2710	CALL message( 'check_parameters', 'PA0129', 1, 2, 0, 6, 0 )
2711	ELSE
2712	DO k = 3, nzt-2
2713	IF ( disturbance_level_t <= zu(k) ) THEN
2714	disturbance_level_ind_t = k
2715	EXIT
2716	ENDIF
2717	ENDDO
2718	ENDIF
2719
2720	!
2721	!-- Check again whether the levels determined this way are ok.
2722	!-- Error may occur at automatic determination and too few grid points in
2723	!-- z-direction.
2724	IF ( disturbance_level_ind_t < disturbance_level_ind_b ) THEN
2725	WRITE( message_string, * ) 'disturbance_level_ind_t = ', &
2726	disturbance_level_ind_t, ' must be >= disturbance_level_b = ', &
2727	disturbance_level_b
2728	CALL message( 'check_parameters', 'PA0130', 1, 2, 0, 6, 0 )
2729	ENDIF
2730
2731	!
2732	!-- Determine the horizontal index range for random perturbations.
2733	!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed
2734	!-- near the inflow and the perturbation area is further limited to ...(1)
2735	!-- after the initial phase of the flow.
2736	dist_nxl = 0; dist_nxr = nx
2737	dist_nys = 0; dist_nyn = ny
2738	IF ( bc_lr /= 'cyclic' ) THEN
2739	IF ( inflow_disturbance_begin == -1 ) THEN
2740	inflow_disturbance_begin = MIN( 10, nx/2 )
2741	ENDIF
2742	IF ( inflow_disturbance_begin < 0 .OR. inflow_disturbance_begin > nx )&
2743	THEN
2744	message_string = 'inflow_disturbance_begin out of range'
2745	CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
2746	ENDIF
2747	IF ( inflow_disturbance_end == -1 ) THEN
2748	inflow_disturbance_end = MIN( 100, 3*nx/4 )
2749	ENDIF
2750	IF ( inflow_disturbance_end < 0 .OR. inflow_disturbance_end > nx ) &
2751	THEN
2752	message_string = 'inflow_disturbance_end out of range'
2753	CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
2754	ENDIF
2755	ELSEIF ( bc_ns /= 'cyclic' ) THEN
2756	IF ( inflow_disturbance_begin == -1 ) THEN
2757	inflow_disturbance_begin = MIN( 10, ny/2 )
2758	ENDIF
2759	IF ( inflow_disturbance_begin < 0 .OR. inflow_disturbance_begin > ny )&
2760	THEN
2761	message_string = 'inflow_disturbance_begin out of range'
2762	CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
2763	ENDIF
2764	IF ( inflow_disturbance_end == -1 ) THEN
2765	inflow_disturbance_end = MIN( 100, 3*ny/4 )
2766	ENDIF
2767	IF ( inflow_disturbance_end < 0 .OR. inflow_disturbance_end > ny ) &
2768	THEN
2769	message_string = 'inflow_disturbance_end out of range'
2770	CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
2771	ENDIF
2772	ENDIF
2773
2774	IF ( bc_lr == 'radiation/dirichlet' ) THEN
2775	dist_nxr = nx - inflow_disturbance_begin
2776	dist_nxl(1) = nx - inflow_disturbance_end
2777	ELSEIF ( bc_lr == 'dirichlet/radiation' ) THEN
2778	dist_nxl = inflow_disturbance_begin
2779	dist_nxr(1) = inflow_disturbance_end
2780	ENDIF
2781	IF ( bc_ns == 'dirichlet/radiation' ) THEN
2782	dist_nyn = ny - inflow_disturbance_begin
2783	dist_nys(1) = ny - inflow_disturbance_end
2784	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
2785	dist_nys = inflow_disturbance_begin
2786	dist_nyn(1) = inflow_disturbance_end
2787	ENDIF
2788
2789	!
2790	!-- A turbulent inflow requires Dirichlet conditions at the respective inflow
2791	!-- boundary (so far, a turbulent inflow is realized from the left side only)
2792	IF ( turbulent_inflow .AND. bc_lr /= 'dirichlet/radiation' ) THEN
2793	message_string = 'turbulent_inflow = .T. requires a Dirichlet ' // &
2794	'condition at the inflow boundary'
2795	CALL message( 'check_parameters', 'PA0133', 1, 2, 0, 6, 0 )
2796	ENDIF
2797
2798	!
2799	!-- In case of turbulent inflow calculate the index of the recycling plane
2800	IF ( turbulent_inflow ) THEN
2801	IF ( recycling_width == 9999999.9 ) THEN
2802	!
2803	!-- Set the default value for the width of the recycling domain
2804	recycling_width = 0.1 * nx * dx
2805	ELSE
2806	IF ( recycling_width < dx .OR. recycling_width > nx * dx ) THEN
2807	WRITE( message_string, * ) 'illegal value for recycling_width:', &
2808	' ', recycling_width
2809	CALL message( 'check_parameters', 'PA0134', 1, 2, 0, 6, 0 )
2810	ENDIF
2811	ENDIF
2812	!
2813	!-- Calculate the index
2814	recycling_plane = recycling_width / dx
2815	ENDIF
2816
2817	!
2818	!-- Check random generator
2819	IF ( random_generator /= 'system-specific' .AND. &
2820	random_generator /= 'numerical-recipes' ) THEN
2821	message_string = 'unknown random generator: random_generator = "' // &
2822	TRIM( random_generator ) // '"'
2823	CALL message( 'check_parameters', 'PA0135', 1, 2, 0, 6, 0 )
2824	ENDIF
2825
2826	!
2827	!-- Determine damping level index for 1D model
2828	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
2829	IF ( damp_level_1d == -1.0 ) THEN
2830	damp_level_1d = zu(nzt+1)
2831	damp_level_ind_1d = nzt + 1
2832	ELSEIF ( damp_level_1d < 0.0 .OR. damp_level_1d > zu(nzt+1) ) THEN
2833	WRITE( message_string, * ) 'damp_level_1d = ', damp_level_1d, &
2834	' must be > 0.0 and < ', zu(nzt+1), '(zu(nzt+1))'
2835	CALL message( 'check_parameters', 'PA0136', 1, 2, 0, 6, 0 )
2836	ELSE
2837	DO k = 1, nzt+1
2838	IF ( damp_level_1d <= zu(k) ) THEN
2839	damp_level_ind_1d = k
2840	EXIT
2841	ENDIF
2842	ENDDO
2843	ENDIF
2844	ENDIF
2845
2846	!
2847	!-- Check some other 1d-model parameters
2848	IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model' .AND. &
2849	TRIM( mixing_length_1d ) /= 'blackadar' ) THEN
2850	message_string = 'mixing_length_1d = "' // TRIM( mixing_length_1d ) // &
2851	'" is unknown'
2852	CALL message( 'check_parameters', 'PA0137', 1, 2, 0, 6, 0 )
2853	ENDIF
2854	IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model' .AND. &
2855	TRIM( dissipation_1d ) /= 'detering' ) THEN
2856	message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) // &
2857	'" is unknown'
2858	CALL message( 'check_parameters', 'PA0138', 1, 2, 0, 6, 0 )
2859	ENDIF
2860
2861	!
2862	!-- Set time for the next user defined restart (time_restart is the
2863	!-- internal parameter for steering restart events)
2864	IF ( restart_time /= 9999999.9 ) THEN
2865	IF ( restart_time > time_since_reference_point ) THEN
2866	time_restart = restart_time
2867	ENDIF
2868	ELSE
2869	!
2870	!-- In case of a restart run, set internal parameter to default (no restart)
2871	!-- if the NAMELIST-parameter restart_time is omitted
2872	time_restart = 9999999.9
2873	ENDIF
2874
2875	!
2876	!-- Set default value of the time needed to terminate a model run
2877	IF ( termination_time_needed == -1.0 ) THEN
2878	IF ( host(1:3) == 'ibm' ) THEN
2879	termination_time_needed = 300.0
2880	ELSE
2881	termination_time_needed = 35.0
2882	ENDIF
2883	ENDIF
2884
2885	!
2886	!-- Check the time needed to terminate a model run
2887	IF ( host(1:3) == 't3e' ) THEN
2888	!
2889	!-- Time needed must be at least 30 seconds on all CRAY machines, because
2890	!-- MPP_TREMAIN gives the remaining CPU time only in steps of 30 seconds
2891	IF ( termination_time_needed <= 30.0 ) THEN
2892	WRITE( message_string, * ) 'termination_time_needed = ', &
2893	termination_time_needed, ' must be > 30.0 on host "', &
2894	TRIM( host ), '"'
2895	CALL message( 'check_parameters', 'PA0139', 1, 2, 0, 6, 0 )
2896	ENDIF
2897	ELSEIF ( host(1:3) == 'ibm' ) THEN
2898	!
2899	!-- On IBM-regatta machines the time should be at least 300 seconds,
2900	!-- because the job time consumed before executing palm (for compiling,
2901	!-- copying of files, etc.) has to be regarded
2902	IF ( termination_time_needed < 300.0 ) THEN
2903	WRITE( message_string, * ) 'termination_time_needed = ', &
2904	termination_time_needed, ' should be >= 300.0 on host "', &
2905	TRIM( host ), '"'
2906	CALL message( 'check_parameters', 'PA0140', 1, 2, 0, 6, 0 )
2907	ENDIF
2908	ENDIF
2909
2910	!
2911	!-- Check pressure gradient conditions
2912	IF ( dp_external .AND. conserve_volume_flow ) THEN
2913	WRITE( message_string, * ) 'Both dp_external and conserve_volume_flow', &
2914	' are .TRUE. but one of them must be .FALSE.'
2915	CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
2916	ENDIF
2917	IF ( dp_external ) THEN
2918	IF ( dp_level_b < zu(nzb) .OR. dp_level_b > zu(nzt) ) THEN
2919	WRITE( message_string, * ) 'dp_level_b = ', dp_level_b, ' is out ', &
2920	' of range'
2921	CALL message( 'check_parameters', 'PA0151', 1, 2, 0, 6, 0 )
2922	ENDIF
2923	IF ( .NOT. ANY( dpdxy /= 0.0 ) ) THEN
2924	WRITE( message_string, * ) 'dp_external is .TRUE. but dpdxy is zero',&
2925	', i.e. the external pressure gradient & will not be applied'
2926	CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
2927	ENDIF
2928	ENDIF
2929	IF ( ANY( dpdxy /= 0.0 ) .AND. .NOT. dp_external ) THEN
2930	WRITE( message_string, * ) 'dpdxy is nonzero but dp_external is ', &
2931	'.FALSE., i.e. the external pressure gradient & will not be applied'
2932	CALL message( 'check_parameters', 'PA0153', 0, 1, 0, 6, 0 )
2933	ENDIF
2934	IF ( conserve_volume_flow ) THEN
2935	IF ( TRIM( conserve_volume_flow_mode ) == 'default' ) THEN
2936	IF ( bc_lr /= 'cyclic' .OR. bc_ns /= 'cyclic' ) THEN
2937	conserve_volume_flow_mode = 'inflow_profile'
2938	ELSE
2939	conserve_volume_flow_mode = 'initial_profiles'
2940	ENDIF
2941	ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles' .AND. &
2942	TRIM( conserve_volume_flow_mode ) /= 'inflow_profile' .AND. &
2943	TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) THEN
2944	WRITE( message_string, * ) 'unknown conserve_volume_flow_mode: ', &
2945	conserve_volume_flow_mode
2946	CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 )
2947	ENDIF
2948	IF ( ( bc_lr /= 'cyclic' .OR. bc_ns /= 'cyclic' ) .AND. &
2949	TRIM( conserve_volume_flow_mode ) /= 'inflow_profile' ) THEN
2950	WRITE( message_string, * ) 'noncyclic boundary conditions ', &
2951	'require & conserve_volume_flow_mode = ''inflow_profile'''
2952	CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 )
2953	ENDIF
2954	IF ( bc_lr == 'cyclic' .AND. bc_ns == 'cyclic' .AND. &
2955	TRIM( conserve_volume_flow_mode ) == 'inflow_profile' ) THEN
2956	WRITE( message_string, * ) 'cyclic boundary conditions ', &
2957	'require & conserve_volume_flow_mode = ''initial_profiles''', &
2958	' or ''bulk_velocity'''
2959	CALL message( 'check_parameters', 'PA0156', 1, 2, 0, 6, 0 )
2960	ENDIF
2961	ENDIF
2962	IF ( ( u_bulk /= 0.0 .OR. v_bulk /= 0.0 ) .AND. &
2963	( .NOT. conserve_volume_flow .OR. &
2964	TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) ) THEN
2965	WRITE( message_string, * ) 'nonzero bulk velocity requires ', &
2966	'conserve_volume_flow = .T. and & ', &
2967	'conserve_volume_flow_mode = ''bulk_velocity'''
2968	CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
2969	ENDIF
2970
2971	!
2972	!-- Check particle attributes
2973	IF ( particle_color /= 'none' ) THEN
2974	IF ( particle_color /= 'absuv' .AND. particle_color /= 'pt*' .AND. &
2975	particle_color /= 'z' ) THEN
2976	message_string = 'illegal value for parameter particle_color: ' // &
2977	TRIM( particle_color)
2978	CALL message( 'check_parameters', 'PA0313', 1, 2, 0, 6, 0 )
2979	ELSE
2980	IF ( color_interval(2) <= color_interval(1) ) THEN
2981	message_string = 'color_interval(2) <= color_interval(1)'
2982	CALL message( 'check_parameters', 'PA0315', 1, 2, 0, 6, 0 )
2983	ENDIF
2984	ENDIF
2985	ENDIF
2986
2987	IF ( particle_dvrpsize /= 'none' ) THEN
2988	IF ( particle_dvrpsize /= 'absw' ) THEN
2989	message_string = 'illegal value for parameter particle_dvrpsize:' // &
2990	' ' // TRIM( particle_color)
2991	CALL message( 'check_parameters', 'PA0314', 1, 2, 0, 6, 0 )
2992	ELSE
2993	IF ( dvrpsize_interval(2) <= dvrpsize_interval(1) ) THEN
2994	message_string = 'dvrpsize_interval(2) <= dvrpsize_interval(1)'
2995	CALL message( 'check_parameters', 'PA0316', 1, 2, 0, 6, 0 )
2996	ENDIF
2997	ENDIF
2998	ENDIF
2999
3000	!
3001	!-- Check &userpar parameters
3002	CALL user_check_parameters
3003
3004
3005	END SUBROUTINE check_parameters

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |