source: palm/tags/release-3.2/SOURCE/check_parameters.f90 @ 258

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

New:
---

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

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

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

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

new inipar-parameter loop_optimization to control the loop optimization method

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

new user interface user_advec_particles

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

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

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

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

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

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

New: wall_fluxes

Changed:


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

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

+age_m in particle_type

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

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

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

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

q is not allowed to become negative (prognostic_equations).

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

d3par-parameter moisture renamed to humidity.

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

vtk directives removed from main program.

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

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

Errors:


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

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

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

Check for possible negative humidities in the initial humidity profile.

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

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

  • Property svn:keywords set to Id
File size: 86.6 KB
Line 
1 SUBROUTINE check_parameters
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: check_parameters.f90 77 2007-03-29 04:26:56Z heinze $
11!
12! 75 2007-03-22 09:54:05Z raasch
13! "by_user" allowed as initializing action, -data_output_ts,
14! leapfrog with non-flat topography not allowed any more, loop_optimization
15! and pt_reference are checked, moisture renamed humidity,
16! output of precipitation amount/rate and roughnes length + check
17! possible negative humidities are avoided in initial profile,
18! dirichlet/neumann changed to dirichlet/radiation, etc.,
19! revision added to run_description_header
20!
21! 20 2007-02-26 00:12:32Z raasch
22! Temperature and humidity gradients at top are now calculated for nzt+1,
23! top_heatflux and respective boundary condition bc_pt_t is checked
24!
25! RCS Log replace by Id keyword, revision history cleaned up
26!
27! Revision 1.61  2006/08/04 14:20:25  raasch
28! do2d_unit and do3d_unit now defined as 2d-arrays, check of
29! use_upstream_for_tke, default value for dt_dopts,
30! generation of file header moved from routines palm and header to here
31!
32! Revision 1.1  1997/08/26 06:29:23  raasch
33! Initial revision
34!
35!
36! Description:
37! ------------
38! Check control parameters and deduce further quantities.
39!------------------------------------------------------------------------------!
40
41    USE arrays_3d
42    USE constants
43    USE control_parameters
44    USE grid_variables
45    USE indices
46    USE model_1d
47    USE netcdf_control
48    USE particle_attributes
49    USE pegrid
50    USE profil_parameter
51    USE statistics
52    USE transpose_indices
53
54    IMPLICIT NONE
55
56    CHARACTER (LEN=1)   ::  sq
57    CHARACTER (LEN=6)   ::  var
58    CHARACTER (LEN=7)   ::  unit
59    CHARACTER (LEN=8)   ::  date
60    CHARACTER (LEN=10)  ::  time
61    CHARACTER (LEN=100) ::  action
62
63    INTEGER ::  i, ilen, intervals, iter, j, k, nnxh, nnyh, position, prec
64    LOGICAL ::  found, ldum
65    REAL    ::  gradient, maxn, maxp
66
67
68!
69!-- Warning, if host is not set
70    IF ( host(1:1) == ' ' )  THEN
71       IF ( myid == 0 )  THEN
72          PRINT*, '+++ WARNING: check_parameters:'
73          PRINT*, '    "host" is not set. Please check that environment', &
74                       ' variable "localhost"'
75          PRINT*, '    is set before running PALM'
76       ENDIF
77    ENDIF
78
79!
80!-- Generate the file header which is used as a header for most of PALM's
81!-- output files
82    CALL DATE_AND_TIME( date, time )
83    run_date = date(7:8)//'-'//date(5:6)//'-'//date(3:4)
84    run_time = time(1:2)//':'//time(3:4)//':'//time(5:6)
85
86    WRITE ( run_description_header, '(A,2X,A,2X,A,A,A,I2.2,2X,A,A,2X,A,1X,A)' )&
87              TRIM( version ), TRIM( revision ), &
88              'run: ', TRIM( run_identifier ), '.', &
89              runnr, 'host: ', TRIM( host ), run_date, run_time
90
91!
92!-- Check the general loop optimization method
93    IF ( loop_optimization == 'default' )  THEN
94       IF ( host(1:3) == 'nec' )  THEN
95          loop_optimization = 'vector'
96       ELSE
97          loop_optimization = 'cache'
98       ENDIF
99    ENDIF
100    IF ( loop_optimization /= 'noopt'  .AND.  loop_optimization /= 'cache' &
101         .AND.  loop_optimization /= 'vector' )  THEN
102       IF ( myid == 0 )  THEN
103          PRINT*, '+++ check_parameters:'
104          PRINT*, '    illegal value given for loop_optimization: ', &
105                  TRIM( loop_optimization )
106       ENDIF
107       CALL local_stop
108    ENDIF
109
110!
111!-- Check topography setting (check for illegal parameter combinations)
112    IF ( topography /= 'flat' )  THEN
113       action = ' '
114       IF ( scalar_advec /= 'pw-scheme' )  THEN
115          WRITE( action, '(A,A)' )  'scalar_advec = ', scalar_advec
116       ENDIF
117       IF ( momentum_advec /= 'pw-scheme' )  THEN
118          WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
119       ENDIF
120       IF ( timestep_scheme(1:8) == 'leapfrog' )  THEN
121          WRITE( action, '(A,A)' )  'timestep_scheme = ', timestep_scheme
122       ENDIF
123       IF ( psolver == 'multigrid'  .OR.  psolver == 'sor' )  THEN
124          WRITE( action, '(A,A)' )  'psolver = ', psolver
125       ENDIF
126       IF ( sloping_surface )  THEN
127          WRITE( action, '(A)' )  'sloping surface = .TRUE.'
128       ENDIF
129       IF ( galilei_transformation )  THEN
130          WRITE( action, '(A)' )  'galilei_transformation = .TRUE.'
131       ENDIF
132       IF ( cloud_physics )  THEN
133          WRITE( action, '(A)' )  'cloud_physics = .TRUE.'
134       ENDIF
135       IF ( cloud_droplets )  THEN
136          WRITE( action, '(A)' )  'cloud_droplets = .TRUE.'
137       ENDIF
138       IF ( humidity )  THEN
139          WRITE( action, '(A)' )  'humidity = .TRUE.'
140       ENDIF
141       IF ( .NOT. prandtl_layer )  THEN
142          WRITE( action, '(A)' )  'prandtl_layer = .FALSE.'
143       ENDIF
144       IF ( action /= ' ' )  THEN
145          IF ( myid == 0 )  THEN
146             PRINT*, '+++ check_parameters:'
147             PRINT*, '    a non-flat topography does not allow ', TRIM( action )
148          ENDIF
149          CALL local_stop
150       ENDIF
151    ENDIF
152!
153!-- Check whether there are any illegal values
154!-- Pressure solver:
155    IF ( psolver /= 'poisfft'  .AND.  psolver /= 'poisfft_hybrid'  .AND. &
156         psolver /= 'sor'  .AND.  psolver /= 'multigrid' )  THEN
157       IF ( myid == 0 )  THEN
158          PRINT*, '+++ check_parameters:'
159          PRINT*, '    unknown solver for perturbation pressure: psolver=', &
160                  psolver
161       ENDIF
162       CALL local_stop
163    ENDIF
164
165#if defined( __parallel )
166    IF ( psolver == 'poisfft_hybrid'  .AND.  pdims(2) /= 1 )  THEN
167       IF ( myid == 0 )  THEN
168          PRINT*, '+++ check_parameters:'
169          PRINT*, '    psolver="', TRIM( psolver ), '" only works for a ', &
170                       '1d domain-decomposition along x'
171          PRINT*, '    please do not set npey/=1 in the parameter file'
172       ENDIF
173       CALL local_stop
174    ENDIF
175    IF ( ( psolver == 'poisfft_hybrid'  .OR.  psolver == 'multigrid' )  .AND.  &
176         ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz ) )  THEN
177       IF ( myid == 0 )  THEN
178          PRINT*, '+++ check_parameters:'
179          PRINT*, '    psolver="', TRIM( psolver ), '" does not work for ', &
180                       'subdomains with unequal size'
181          PRINT*, '    please set grid_matching = ''strict'' in the parameter',&
182                       ' file'
183       ENDIF
184       CALL local_stop
185    ENDIF
186#else
187    IF ( psolver == 'poisfft_hybrid' )  THEN
188       IF ( myid == 0 )  THEN
189          PRINT*, '+++ check_parameters:'
190          PRINT*, '    psolver="', TRIM( psolver ), '" only works for a ', &
191                       'parallel environment'
192       ENDIF
193       CALL local_stop
194    ENDIF
195#endif
196
197    IF ( psolver == 'multigrid' )  THEN
198       IF ( cycle_mg == 'w' )  THEN
199          gamma_mg = 2
200       ELSEIF ( cycle_mg == 'v' )  THEN
201          gamma_mg = 1
202       ELSE
203          IF ( myid == 0 )  THEN
204             PRINT*, '+++ check_parameters:'
205             PRINT*, '    unknown multigrid cycle: cycle_mg=', cycle_mg
206          ENDIF
207          CALL local_stop
208       ENDIF
209    ENDIF
210
211    IF ( fft_method /= 'singleton-algorithm'  .AND.  &
212         fft_method /= 'temperton-algorithm'  .AND.  &
213         fft_method /= 'system-specific' )  THEN
214       IF ( myid == 0 )  THEN
215          PRINT*, '+++ check_parameters:'
216          PRINT*, '    unknown fft-algorithm: fft_method=', fft_method
217       ENDIF
218       CALL local_stop
219    ENDIF
220
221!
222!-- Advection schemes:
223    IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ups-scheme' ) &
224    THEN
225       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  unknown advection ', &
226                                 'scheme: momentum_advec=', momentum_advec
227       CALL local_stop
228    ENDIF
229    IF ( ( momentum_advec == 'ups-scheme'  .OR.  scalar_advec == 'ups-scheme' )&
230                                      .AND.  timestep_scheme /= 'euler' )  THEN
231       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  momentum_advec=', &
232                                 momentum_advec, ' is not allowed with ', &
233                                 'timestep_scheme=', timestep_scheme
234       CALL local_stop
235    ENDIF
236
237    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'bc-scheme'  .AND.&
238         scalar_advec /= 'ups-scheme' )  THEN
239       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  unknown advection ', &
240                                 'scheme: scalar_advec=', scalar_advec
241       CALL local_stop
242    ENDIF
243
244    IF ( use_sgs_for_particles  .AND.  .NOT. use_upstream_for_tke )  THEN
245       use_upstream_for_tke = .TRUE.
246       IF ( myid == 0 )  THEN
247          PRINT*, '+++ WARNING: check_parameters:  use_upstream_for_tke set ', &
248                       '.TRUE. because use_sgs_for_particles = .TRUE.'
249       ENDIF
250    ENDIF
251
252    IF ( use_upstream_for_tke  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
253       IF ( myid == 0 )  THEN
254          PRINT*, '+++ check_parameters:  use_upstream_for_tke = .TRUE. ', &
255                       'not allowed with timestep_scheme=', timestep_scheme
256       ENDIF
257       CALL local_stop
258    ENDIF
259
260!
261!-- Timestep schemes:
262    SELECT CASE ( TRIM( timestep_scheme ) )
263
264       CASE ( 'euler' )
265          intermediate_timestep_count_max = 1
266          asselin_filter_factor           = 0.0
267
268       CASE ( 'leapfrog', 'leapfrog+euler' )
269          intermediate_timestep_count_max = 1
270
271       CASE ( 'runge-kutta-2' )
272          intermediate_timestep_count_max = 2
273          asselin_filter_factor           = 0.0
274
275       CASE ( 'runge-kutta-3' )
276          intermediate_timestep_count_max = 3
277          asselin_filter_factor           = 0.0
278
279       CASE DEFAULT
280          IF ( myid == 0 )  PRINT*, '+++ check_parameters:  unknown timestep ',&
281                                    'scheme: timestep_scheme=', timestep_scheme
282          CALL local_stop
283
284    END SELECT
285
286    IF ( scalar_advec == 'ups-scheme'  .AND.  timestep_scheme(1:5) == 'runge' )&
287    THEN
288       IF ( myid == 0 )  THEN
289          PRINT*, '+++ check_parameters:  scalar advection scheme "', &
290                                          TRIM( scalar_advec ), '"'
291          PRINT*, '    does not work with timestep_scheme "', &
292                                          TRIM( timestep_scheme ), '"'
293       ENDIF
294       CALL local_stop
295    ENDIF
296
297    IF ( momentum_advec /= 'pw-scheme' .AND. timestep_scheme(1:5) == 'runge' ) &
298    THEN
299       IF ( myid == 0 )  THEN
300          PRINT*, '+++ check_parameters:  momentum advection scheme "', &
301                                          TRIM( momentum_advec ), '"'
302          PRINT*, '    does not work with timestep_scheme "', &
303                                          TRIM( timestep_scheme ), '"'
304       ENDIF
305       CALL local_stop
306    ENDIF
307
308
309    IF ( initializing_actions == ' ' )  THEN
310       IF ( myid == 0 )  THEN
311          PRINT*, '+++ check parameters:'
312          PRINT*, '    no value found for initializing_actions'
313       ENDIF
314       CALL local_stop
315    ENDIF
316
317    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
318!
319!--    No model continuation run; several initialising actions are possible
320       action = initializing_actions
321       DO WHILE ( TRIM( action ) /= '' )
322          position = INDEX( action, ' ' )
323          SELECT CASE ( action(1:position-1) )
324
325             CASE ( 'set_constant_profiles', 'set_1d-model_profiles', &
326                    'by_user', 'initialize_vortex',     'initialize_ptanom' )
327                action = action(position+1:)
328
329             CASE DEFAULT
330                IF ( myid == 0 )  PRINT*, '+++ check_parameters: initializi', &
331                               'ng_action unkown or not allowed: action = "', &
332                               TRIM(action), '"'
333                CALL local_stop
334
335          END SELECT
336       ENDDO
337    ENDIF
338    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND. &
339         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
340       IF ( myid == 0 )  PRINT*, '+++ check_parameters: initializing_actions', &
341          '"set_constant_profiles" and "set_1d-model_profiles" are not', &
342          ' allowed simultaneously'
343       CALL local_stop
344    ENDIF
345    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND. &
346         INDEX( initializing_actions, 'by_user' ) /= 0 )  THEN
347       IF ( myid == 0 )  PRINT*, '+++ check_parameters: initializing_actions', &
348          '"set_constant_profiles" and "by_user" are not', &
349          ' allowed simultaneously'
350       CALL local_stop
351    ENDIF
352    IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND. &
353         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
354       IF ( myid == 0 )  PRINT*, '+++ check_parameters: initializing_actions', &
355          '"by_user" and "set_1d-model_profiles" are not', &
356          ' allowed simultaneously'
357       CALL local_stop
358    ENDIF
359
360    IF ( cloud_physics  .AND.  .NOT. humidity )  THEN
361       IF ( myid == 0 )  PRINT*, '+++ check_parameters: cloud_physics =', &
362                                 cloud_physics, ' is not allowed with ',  &
363                                 'humidity =', humidity
364       CALL local_stop
365    ENDIF
366
367    IF ( precipitation  .AND.  .NOT.  cloud_physics )  THEN
368       IF ( myid == 0 )  PRINT*, '+++ check_parameters: precipitation =', &
369                                 precipitation, ' is not allowed with ',  &
370                                 'cloud_physics =', cloud_physics
371       CALL local_stop
372    ENDIF
373
374    IF ( humidity  .AND.  sloping_surface )  THEN
375       IF ( myid == 0 )  PRINT*, '+++ check_parameters: humidity = TRUE', &
376                                 'and hang = TRUE are not',               &
377                                 ' allowed simultaneously' 
378       CALL local_stop       
379    ENDIF
380
381    IF ( humidity  .AND.  scalar_advec == 'ups-scheme' )  THEN
382       IF ( myid == 0 )  PRINT*, '+++ check_parameters: UPS-scheme', &
383                                 'is not implemented for humidity'
384       CALL local_stop       
385    ENDIF
386
387    IF ( passive_scalar  .AND.  humidity )  THEN
388       IF ( myid == 0 )  PRINT*, '+++ check_parameters: humidity = TRUE and', &
389                                 'passive_scalar = TRUE is not allowed ',     &
390                                 'simultaneously'
391       CALL local_stop
392    ENDIF
393
394    IF ( passive_scalar  .AND.  scalar_advec == 'ups-scheme' )  THEN
395       IF ( myid == 0 )  PRINT*, '+++ check_parameters: UPS-scheme', &
396                                 'is not implemented for passive_scalar'
397       CALL local_stop       
398    ENDIF
399
400    IF ( grid_matching /= 'strict'  .AND.  grid_matching /= 'match' )  THEN
401       IF ( myid == 0 )  PRINT*, '+++ check_parameters: illegal value "', &
402                                 TRIM( grid_matching ),                   &
403                                 '" found for parameter grid_matching'
404       CALL local_stop       
405    ENDIF
406
407!
408!-- In case of no model continuation run, check initialising parameters and
409!-- deduce further quantities
410    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
411
412!
413!--    Initial profiles for 1D and 3D model, respectively
414       u_init  = ug_surface
415       v_init  = vg_surface
416       pt_init = pt_surface
417       IF ( humidity )        q_init = q_surface
418       IF ( passive_scalar )  q_init = s_surface
419
420!
421!--
422!--    If required, compute initial profile of the geostrophic wind
423!--    (component ug)
424       i = 1
425       gradient = 0.0
426       ug_vertical_gradient_level_ind(1) = 0
427       ug(0) = ug_surface
428       DO  k = 1, nzt+1
429          IF ( ug_vertical_gradient_level(i) < zu(k)  .AND. &
430               ug_vertical_gradient_level(i) >= 0.0 )  THEN
431             gradient = ug_vertical_gradient(i) / 100.0
432             ug_vertical_gradient_level_ind(i) = k - 1
433             i = i + 1
434             IF ( i > 10 )  THEN
435                IF ( myid == 0 )  THEN
436                   PRINT*, '+++ check_parameters: upper bound 10 of array', &
437                           ' "ug_vertical_gradient_level_ind" exceeded'
438                ENDIF
439                CALL local_stop
440             ENDIF
441          ENDIF
442          IF ( gradient /= 0.0 )  THEN
443             IF ( k /= 1 )  THEN
444                ug(k) = ug(k-1) + dzu(k) * gradient
445             ELSE
446                ug(k) = ug_surface + 0.5 * dzu(k) * gradient
447             ENDIF
448          ELSE
449             ug(k) = ug(k-1)
450          ENDIF
451       ENDDO
452
453       u_init = ug
454
455!
456!--    In case of no given gradients for ug, choose a vanishing gradient
457       IF ( ug_vertical_gradient_level(1) == -1.0 )  THEN
458          ug_vertical_gradient_level(1) = 0.0
459       ENDIF 
460
461!
462!--
463!--    If required, compute initial profile of the geostrophic wind
464!--    (component vg)
465       i = 1
466       gradient = 0.0
467       vg_vertical_gradient_level_ind(1) = 0
468       vg(0) = vg_surface
469       DO  k = 1, nzt+1
470          IF ( vg_vertical_gradient_level(i) < zu(k)  .AND. &
471               vg_vertical_gradient_level(i) >= 0.0 )  THEN
472             gradient = vg_vertical_gradient(i) / 100.0
473             vg_vertical_gradient_level_ind(i) = k - 1
474             i = i + 1
475             IF ( i > 10 )  THEN
476                IF ( myid == 0 )  THEN
477                   PRINT*, '+++ check_parameters: upper bound 10 of array', &
478                           ' "vg_vertical_gradient_level_ind" exceeded'
479                ENDIF
480                CALL local_stop
481             ENDIF
482          ENDIF
483          IF ( gradient /= 0.0 )  THEN
484             IF ( k /= 1 )  THEN
485                vg(k) = vg(k-1) + dzu(k) * gradient
486             ELSE
487                vg(k) = vg_surface + 0.5 * dzu(k) * gradient
488             ENDIF
489          ELSE
490             vg(k) = vg(k-1)
491          ENDIF
492       ENDDO
493
494       v_init = vg
495 
496!
497!--    In case of no given gradients for vg, choose a vanishing gradient
498       IF ( vg_vertical_gradient_level(1) == -1.0 )  THEN
499          vg_vertical_gradient_level(1) = 0.0
500       ENDIF
501
502!
503!--    If required, compute initial temperature profile using the given
504!--    temperature gradient
505       i = 1
506       gradient = 0.0
507       pt_vertical_gradient_level_ind(1) = 0
508       DO  k = 1, nzt+1
509          IF ( pt_vertical_gradient_level(i) < zu(k)  .AND. &
510               pt_vertical_gradient_level(i) >= 0.0 )  THEN
511             gradient = pt_vertical_gradient(i) / 100.0
512             pt_vertical_gradient_level_ind(i) = k - 1
513             i = i + 1
514             IF ( i > 10 )  THEN
515                IF ( myid == 0 )  THEN
516                   PRINT*, '+++ check_parameters: upper bound 10 of array', &
517                           ' "pt_vertical_gradient_level_ind" exceeded'
518                ENDIF
519                CALL local_stop
520             ENDIF
521          ENDIF
522          IF ( gradient /= 0.0 )  THEN
523             IF ( k /= 1 )  THEN
524                pt_init(k) = pt_init(k-1) + dzu(k) * gradient
525             ELSE
526                pt_init(k) = pt_init(k-1) + 0.5 * dzu(k) * gradient
527             ENDIF
528          ELSE
529             pt_init(k) = pt_init(k-1)
530          ENDIF
531       ENDDO
532
533!
534!--    In case of no given temperature gradients, choose gradient of neutral
535!--    stratification
536       IF ( pt_vertical_gradient_level(1) == -1.0 )  THEN
537          pt_vertical_gradient_level(1) = 0.0
538       ENDIF
539
540!
541!--    Store temperature gradient at the top boundary for possile Neumann
542!--    boundary condition
543       bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
544
545!
546!--    If required, compute initial humidity or scalar profile using the given
547!--    humidity/scalar gradient. In case of scalar transport, initially store
548!--    values of the scalar parameters on humidity parameters
549       IF ( passive_scalar )  THEN
550          bc_q_b                    = bc_s_b
551          bc_q_t                    = bc_s_t
552          q_surface                 = s_surface
553          q_surface_initial_change  = s_surface_initial_change
554          q_vertical_gradient       = s_vertical_gradient
555          q_vertical_gradient_level = s_vertical_gradient_level
556          surface_waterflux         = surface_scalarflux
557       ENDIF
558
559       IF ( humidity  .OR.  passive_scalar )  THEN
560
561          i = 1
562          gradient = 0.0
563          q_vertical_gradient_level_ind(1) = 0
564          DO  k = 1, nzt+1
565             IF ( q_vertical_gradient_level(i) < zu(k)  .AND. &
566                  q_vertical_gradient_level(i) >= 0.0 )  THEN
567                gradient = q_vertical_gradient(i) / 100.0
568                q_vertical_gradient_level_ind(i) = k - 1
569                i = i + 1
570                IF ( i > 10 )  THEN
571                   IF ( myid == 0 )  THEN
572                      PRINT*, '+++ check_parameters: upper bound 10 of arr', &
573                              'ay "q_vertical_gradient_level_ind" exceeded'
574                   ENDIF
575                   CALL local_stop
576                ENDIF
577             ENDIF
578             IF ( gradient /= 0.0 )  THEN
579                IF ( k /= 1 )  THEN
580                   q_init(k) = q_init(k-1) + dzu(k) * gradient
581                ELSE
582                   q_init(k) = q_init(k-1) + 0.5 * dzu(k) * gradient
583                ENDIF
584             ELSE
585                q_init(k) = q_init(k-1)
586             ENDIF
587!
588!--          Avoid negative humidities
589             IF ( q_init(k) < 0.0 )  THEN
590                q_init(k) = 0.0
591             ENDIF
592          ENDDO
593
594!
595!--       In case of no given humidity gradients, choose zero gradient
596!--       conditions
597          IF ( q_vertical_gradient_level(1) == -1.0 )  THEN
598             q_vertical_gradient_level(1) = 0.0
599          ENDIF
600
601!
602!--       Store humidity gradient at the top boundary for possile Neumann
603!--       boundary condition
604          bc_q_t_val = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
605
606       ENDIF
607
608    ENDIF
609
610!
611!-- Compute Coriolis parameter
612    f  = 2.0 * omega * SIN( phi / 180.0 * pi )
613    fs = 2.0 * omega * COS( phi / 180.0 * pi )
614
615!
616!-- Reference temperature to be used in buoyancy terms
617    IF ( pt_reference /= 9999999.9 )  use_pt_reference = .TRUE.
618
619!
620!-- In case of a given slope, compute the relevant quantities
621    IF ( alpha_surface /= 0.0 )  THEN
622       IF ( ABS( alpha_surface ) > 90.0 )  THEN
623          IF ( myid == 0 )  PRINT*, '+++ check_parameters: ABS( alpha_surface',&
624                                    '=', alpha_surface, ' ) must be < 90.0'
625          CALL local_stop
626       ENDIF
627       sloping_surface = .TRUE.
628       cos_alpha_surface = COS( alpha_surface / 180.0 * pi )
629       sin_alpha_surface = SIN( alpha_surface / 180.0 * pi )
630    ENDIF
631
632!
633!-- Check time step and cfl_factor
634    IF ( dt /= -1.0 )  THEN
635       IF ( dt <= 0.0  .AND.  dt /= -1.0 )  THEN
636          IF ( myid == 0 )  PRINT*, '+++ check_parameters:  dt=', dt, ' <= 0.0'
637          CALL local_stop
638       ENDIF
639       dt_3d = dt
640       dt_fixed = .TRUE.
641    ENDIF
642
643    IF ( cfl_factor <= 0.0  .OR.  cfl_factor > 1.0 )  THEN
644       IF ( cfl_factor == -1.0 )  THEN
645          IF ( momentum_advec == 'ups-scheme'  .OR.  &
646               scalar_advec == 'ups-scheme' )  THEN
647             cfl_factor = 0.8
648          ELSE
649             IF ( timestep_scheme == 'runge-kutta-2' )  THEN
650                cfl_factor = 0.8
651             ELSEIF ( timestep_scheme == 'runge-kutta-3' )  THEN
652                cfl_factor = 0.9
653             ELSE
654                cfl_factor = 0.1
655             ENDIF
656          ENDIF
657       ELSE
658          IF ( myid == 0 )  THEN
659             PRINT*, '+++ check_parameters: cfl_factor=', cfl_factor, &
660                         ' out of range'
661             PRINT*, '+++                   0.0 < cfl_factor <= 1.0 is required'
662          ENDIF
663          CALL local_stop
664       ENDIF
665    ENDIF
666
667!
668!-- Store simulated time at begin
669    simulated_time_at_begin = simulated_time
670
671!
672!-- Set wind speed in the Galilei-transformed system
673    IF ( galilei_transformation )  THEN
674       IF ( use_ug_for_galilei_tr .AND.                &
675            ug_vertical_gradient_level(1) == 0.0 .AND. & 
676            vg_vertical_gradient_level(1) == 0.0 )  THEN
677          u_gtrans = ug_surface
678          v_gtrans = vg_surface
679       ELSEIF ( use_ug_for_galilei_tr .AND.                &
680                ug_vertical_gradient_level(1) /= 0.0 )  THEN
681          IF ( myid == 0 )  THEN
682             PRINT*, '+++ check_parameters:'
683             PRINT*, '    baroclinicity (ug) not allowed'
684             PRINT*, '    simultaneously with galilei transformation'
685          ENDIF
686          CALL local_stop
687       ELSEIF ( use_ug_for_galilei_tr .AND.                &
688                vg_vertical_gradient_level(1) /= 0.0 )  THEN
689          IF ( myid == 0 )  THEN
690             PRINT*, '+++ check_parameters:'
691             PRINT*, '    baroclinicity (vg) not allowed'
692             PRINT*, '    simultaneously with galilei transformation'
693          ENDIF
694          CALL local_stop
695       ELSE
696          IF ( myid == 0 )  THEN
697             PRINT*, '+++ WARNING: check_parameters:'
698             PRINT*, '    variable translation speed used for galilei-tran' // &
699                          'sformation, which'
700             PRINT*, '    may cause instabilities in stably stratified regions'
701          ENDIF
702       ENDIF
703    ENDIF
704
705!
706!-- In case of using a prandtl-layer, calculated (or prescribed) surface
707!-- fluxes have to be used in the diffusion-terms
708    IF ( prandtl_layer )  use_surface_fluxes = .TRUE.
709
710!
711!-- Check boundary conditions and set internal variables:
712!-- Lateral boundary conditions
713    IF ( bc_lr /= 'cyclic'  .AND.  bc_lr /= 'dirichlet/radiation'  .AND. &
714         bc_lr /= 'radiation/dirichlet' )  THEN
715       IF ( myid == 0 )  THEN
716          PRINT*, '+++ check_parameters:'
717          PRINT*, '    unknown boundary condition: bc_lr = ', bc_lr
718       ENDIF
719       CALL local_stop
720    ENDIF
721    IF ( bc_ns /= 'cyclic'  .AND.  bc_ns /= 'dirichlet/radiation'  .AND. &
722         bc_ns /= 'radiation/dirichlet' )  THEN
723       IF ( myid == 0 )  THEN
724          PRINT*, '+++ check_parameters:'
725          PRINT*, '    unknown boundary condition: bc_ns = ', bc_ns
726       ENDIF
727       CALL local_stop
728    ENDIF
729
730!
731!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
732!-- Willimas advection scheme. Several schemes and tools do not work with
733!-- non-cyclic boundary conditions.
734    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
735       IF ( psolver /= 'multigrid' )  THEN
736          IF ( myid == 0 )  THEN
737             PRINT*, '+++ check_parameters:'
738             PRINT*, '    non-cyclic lateral boundaries do not allow', &
739                          ' psolver = ', psolver 
740          ENDIF
741          CALL local_stop
742       ENDIF
743       IF ( momentum_advec /= 'pw-scheme' )  THEN
744          IF ( myid == 0 )  THEN
745             PRINT*, '+++ check_parameters:'
746             PRINT*, '    non-cyclic lateral boundaries do not allow', &
747                          ' momentum_advec = ', momentum_advec 
748          ENDIF
749          CALL local_stop
750       ENDIF
751       IF ( scalar_advec /= 'pw-scheme' )  THEN
752          IF ( myid == 0 )  THEN
753             PRINT*, '+++ check_parameters:'
754             PRINT*, '    non-cyclic lateral boundaries do not allow', &
755                          ' scalar_advec = ', scalar_advec 
756          ENDIF
757          CALL local_stop
758       ENDIF
759       IF ( galilei_transformation )  THEN
760          IF ( myid == 0 )  THEN
761             PRINT*, '+++ check_parameters:'
762             PRINT*, '    non-cyclic lateral boundaries do not allow', &
763                          ' galilei_transformation = .T.' 
764          ENDIF
765          CALL local_stop
766       ENDIF
767!       IF ( conserve_volume_flow )  THEN
768!          IF ( myid == 0 )  THEN
769!             PRINT*, '+++ check_parameters:'
770!             PRINT*, '    non-cyclic lateral boundaries do not allow', &
771!                          ' conserve_volume_flow = .T.'
772!          ENDIF
773!          CALL local_stop
774!       ENDIF
775    ENDIF
776
777!
778!-- Bottom boundary condition for the turbulent Kinetic energy
779    IF ( bc_e_b == 'neumann' )  THEN
780       ibc_e_b = 1
781       IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
782          IF ( myid == 0 )  THEN
783             PRINT*, '+++ WARNING: check_parameters:'
784             PRINT*, '    adjust_mixing_length = TRUE and bc_e_b = ', bc_e_b
785          ENDIF
786       ENDIF
787    ELSEIF ( bc_e_b == '(u*)**2+neumann' )  THEN
788       ibc_e_b = 2
789       IF ( .NOT. adjust_mixing_length  .AND.  prandtl_layer )  THEN
790          IF ( myid == 0 )  THEN
791             PRINT*, '+++ WARNING: check_parameters:'
792             PRINT*, '    adjust_mixing_length = FALSE and bc_e_b = ', bc_e_b
793          ENDIF
794       ENDIF
795       IF ( .NOT. prandtl_layer )  THEN
796          bc_e_b = 'neumann'
797          ibc_e_b = 1
798          IF ( myid == 0 )  THEN
799             PRINT*, '+++ WARNING: check_parameters:'
800             PRINT*, '    boundary condition bc_e_b changed to "', bc_e_b, '"'
801          ENDIF
802       ENDIF
803    ELSE
804       IF ( myid == 0 )  THEN
805          PRINT*, '+++ check_parameters:'
806          PRINT*, '    unknown boundary condition: bc_e_b = ', bc_e_b
807       ENDIF
808       CALL local_stop
809    ENDIF
810
811!
812!-- Boundary conditions for perturbation pressure
813    IF ( bc_p_b == 'dirichlet' )  THEN
814       ibc_p_b = 0
815    ELSEIF ( bc_p_b == 'neumann' )  THEN
816       ibc_p_b = 1
817    ELSEIF ( bc_p_b == 'neumann+inhomo' )  THEN
818       ibc_p_b = 2
819    ELSE
820       IF ( myid == 0 )  THEN
821          PRINT*, '+++ check_parameters:'
822          PRINT*, '    unknown boundary condition: bc_p_b = ', bc_p_b
823       ENDIF
824       CALL local_stop
825    ENDIF
826    IF ( ibc_p_b == 2  .AND.  .NOT. prandtl_layer )  THEN
827       IF ( myid == 0 )  THEN
828          PRINT*, '+++ check_parameters:'
829          PRINT*, '    boundary condition: bc_p_b = ', TRIM( bc_p_b ), &
830                       ' not allowed with'
831          PRINT*, '    prandtl_layer = .FALSE.' 
832       ENDIF
833       CALL local_stop
834    ENDIF
835    IF ( bc_p_t == 'dirichlet' )  THEN
836       ibc_p_t = 0
837    ELSEIF ( bc_p_t == 'neumann' )  THEN
838       ibc_p_t = 1
839    ELSE
840       IF ( myid == 0 )  THEN
841          PRINT*, '+++ check_parameters:'
842          PRINT*, '    unknown boundary condition: bc_p_t = ', bc_p_t
843       ENDIF
844       CALL local_stop
845    ENDIF
846
847!
848!-- Boundary conditions for potential temperature
849    IF ( bc_pt_b == 'dirichlet' )  THEN
850       ibc_pt_b = 0
851    ELSEIF ( bc_pt_b == 'neumann' )  THEN
852       ibc_pt_b = 1
853    ELSE
854       IF ( myid == 0 )  THEN
855          PRINT*, '+++ check_parameters:'
856          PRINT*, '    unknown boundary condition: bc_pt_b = ', bc_pt_b
857       ENDIF
858       CALL local_stop
859    ENDIF
860    IF ( bc_pt_t == 'dirichlet' )  THEN
861       ibc_pt_t = 0
862    ELSEIF ( bc_pt_t == 'neumann' )  THEN
863       ibc_pt_t = 1
864    ELSEIF ( bc_pt_t == 'initial_gradient' )  THEN
865       ibc_pt_t = 2
866    ELSE
867       IF ( myid == 0 )  THEN
868          PRINT*, '+++ check_parameters:'
869          PRINT*, '    unknown boundary condition: bc_pt_t = ', bc_pt_t
870       ENDIF
871       CALL local_stop
872    ENDIF
873
874    IF ( surface_heatflux == 9999999.9 )  constant_heatflux     = .FALSE.
875    IF ( top_heatflux     == 9999999.9 )  constant_top_heatflux = .FALSE.
876
877!
878!-- A given surface temperature implies Dirichlet boundary condition for
879!-- temperature. In this case specification of a constant heat flux is
880!-- forbidden.
881    IF ( ibc_pt_b == 0  .AND.   constant_heatflux  .AND. &
882         surface_heatflux /= 0.0 )  THEN
883       IF ( myid == 0 )  THEN
884          PRINT*, '+++ check_parameters:'
885          PRINT*, '    boundary_condition: bc_pt_b = ', bc_pt_b
886          PRINT*, '    is not allowed with constant_heatflux = .TRUE.'
887       ENDIF
888       CALL local_stop
889    ENDIF
890    IF ( constant_heatflux  .AND.  pt_surface_initial_change /= 0.0 )  THEN
891       IF ( myid == 0 )  THEN
892          PRINT*, '+++ check_parameters: constant_heatflux = .TRUE. is not'
893          PRINT*, '    allowed with pt_surface_initial_change (/=0) = ', &
894                  pt_surface_initial_change
895       ENDIF
896       CALL local_stop
897    ENDIF
898
899!
900!-- A given temperature at the top implies Dirichlet boundary condition for
901!-- temperature. In this case specification of a constant heat flux is
902!-- forbidden.
903    IF ( ibc_pt_t == 0  .AND.   constant_top_heatflux  .AND. &
904         top_heatflux /= 0.0 )  THEN
905       IF ( myid == 0 )  THEN
906          PRINT*, '+++ check_parameters:'
907          PRINT*, '    boundary_condition: bc_pt_t = ', bc_pt_t
908          PRINT*, '    is not allowed with constant_top_heatflux = .TRUE.'
909       ENDIF
910       CALL local_stop
911    ENDIF
912
913!
914!-- In case of humidity or passive scalar, set boundary conditions for total
915!-- water content / scalar
916    IF ( humidity  .OR.  passive_scalar ) THEN
917       IF ( humidity )  THEN
918          sq = 'q'
919       ELSE
920          sq = 's'
921       ENDIF
922       IF ( bc_q_b == 'dirichlet' )  THEN
923          ibc_q_b = 0
924       ELSEIF ( bc_q_b == 'neumann' )  THEN
925          ibc_q_b = 1
926       ELSE
927          IF ( myid == 0 )  THEN
928             PRINT*, '+++ check_parameters:'
929             PRINT*, '    unknown boundary condition: bc_', sq, '_b = ', bc_q_b
930          ENDIF
931          CALL local_stop
932       ENDIF
933       IF ( bc_q_t == 'dirichlet' )  THEN
934          ibc_q_t = 0
935       ELSEIF ( bc_q_t == 'neumann' )  THEN
936          ibc_q_t = 1
937       ELSE
938          IF ( myid == 0 )  THEN
939             PRINT*, '+++ check_parameters:'
940             PRINT*, '    unknown boundary condition: bc_', sq, '_t = ', bc_q_t
941          ENDIF
942          CALL local_stop
943       ENDIF
944
945       IF ( surface_waterflux == 0.0 )  constant_waterflux = .FALSE.
946
947!
948!--    A given surface humidity implies Dirichlet boundary condition for
949!--    humidity. In this case specification of a constant water flux is
950!--    forbidden.
951       IF ( ibc_q_b == 0  .AND.  constant_waterflux )  THEN
952          IF ( myid == 0 )  THEN
953             PRINT*, '+++ check_parameters:'
954             PRINT*, '    boundary_condition: bc_', sq, '_b = ', bc_q_b
955             PRINT*, '    is not allowed with prescribed surface flux'
956          ENDIF
957          CALL local_stop
958       ENDIF
959       IF ( constant_waterflux  .AND.  q_surface_initial_change /= 0.0 )  THEN
960          IF ( myid == 0 )  THEN
961             PRINT*, '+++ check_parameters: a prescribed surface flux is not'
962             PRINT*, '    allowed with ', sq, '_surface_initial_change (/=0)', &
963                     ' = ', q_surface_initial_change
964          ENDIF
965          CALL local_stop
966       ENDIF
967       
968    ENDIF
969
970!
971!-- Boundary conditions for horizontal components of wind speed
972    IF ( bc_uv_b == 'dirichlet' )  THEN
973       ibc_uv_b = 0
974    ELSEIF ( bc_uv_b == 'neumann' )  THEN
975       ibc_uv_b = 1
976       IF ( prandtl_layer )  THEN
977          IF ( myid == 0 )  THEN
978             PRINT*, '+++ check_parameters:'
979             PRINT*, '    boundary condition: bc_uv_b = ', TRIM( bc_uv_b ), &
980                          ' is not allowed with'
981             PRINT*, '    prandtl_layer = .TRUE.' 
982          ENDIF
983          CALL local_stop
984       ENDIF
985    ELSE
986       IF ( myid == 0 )  THEN
987          PRINT*, '+++ check_parameters:'
988          PRINT*, '    unknown boundary condition: bc_uv_b = ', bc_uv_b
989       ENDIF
990       CALL local_stop
991    ENDIF
992    IF ( bc_uv_t == 'dirichlet' )  THEN
993       ibc_uv_t = 0
994    ELSEIF ( bc_uv_t == 'neumann' )  THEN
995       ibc_uv_t = 1
996    ELSE
997       IF ( myid == 0 )  THEN
998          PRINT*, '+++ check_parameters:'
999          PRINT*, '    unknown boundary condition: bc_uv_t = ', bc_uv_t
1000       ENDIF
1001       CALL local_stop
1002    ENDIF
1003
1004!
1005!-- Compute and check, respectively, the Rayleigh Damping parameter
1006    IF ( rayleigh_damping_factor == -1.0 )  THEN
1007       IF ( momentum_advec == 'ups-scheme' )  THEN
1008          rayleigh_damping_factor = 0.01
1009       ELSE
1010          rayleigh_damping_factor = 0.0
1011       ENDIF
1012    ELSE
1013       IF ( rayleigh_damping_factor < 0.0 .OR. rayleigh_damping_factor > 1.0 ) &
1014       THEN
1015          IF ( myid == 0 )  THEN
1016             PRINT*, '+++ check_parameters:'
1017             PRINT*, '    rayleigh_damping_factor = ', rayleigh_damping_factor,&
1018                          ' out of range [0.0,1.0]'
1019          ENDIF
1020          CALL local_stop
1021       ENDIF
1022    ENDIF
1023
1024    IF ( rayleigh_damping_height == -1.0 )  THEN
1025       rayleigh_damping_height = 0.66666666666 * zu(nzt)
1026    ELSE
1027       IF ( rayleigh_damping_height < 0.0  .OR. &
1028            rayleigh_damping_height > zu(nzt) )  THEN
1029          IF ( myid == 0 )  THEN
1030             PRINT*, '+++ check_parameters:'
1031             PRINT*, '    rayleigh_damping_height = ', rayleigh_damping_height,&
1032                          ' out of range [0.0,', zu(nzt), ']'
1033          ENDIF
1034          CALL local_stop
1035       ENDIF
1036    ENDIF
1037
1038!
1039!-- Check limiters for Upstream-Spline scheme
1040    IF ( overshoot_limit_u < 0.0  .OR.  overshoot_limit_v < 0.0  .OR.  &
1041         overshoot_limit_w < 0.0  .OR.  overshoot_limit_pt < 0.0  .OR. &
1042         overshoot_limit_e < 0.0 )  THEN
1043       IF ( myid == 0 )  THEN
1044          PRINT*, '+++ check_parameters:'
1045          PRINT*, '    overshoot_limit_... < 0.0 is not allowed'
1046       ENDIF
1047       CALL local_stop
1048    ENDIF
1049    IF ( ups_limit_u < 0.0 .OR. ups_limit_v < 0.0 .OR. ups_limit_w < 0.0 .OR. &
1050         ups_limit_pt < 0.0 .OR. ups_limit_e < 0.0 )  THEN
1051       IF ( myid == 0 )  THEN
1052          PRINT*, '+++ check_parameters:'
1053          PRINT*, '    ups_limit_... < 0.0 is not allowed'
1054       ENDIF
1055       CALL local_stop
1056    ENDIF
1057
1058!
1059!-- Check number of chosen statistic regions. More than 10 regions are not
1060!-- allowed, because so far no more than 10 corresponding output files can
1061!-- be opened (cf. check_open)
1062    IF ( statistic_regions > 9  .OR.  statistic_regions < 0 )  THEN
1063       IF ( myid == 0 )  THEN
1064          PRINT*, '+++ check_parameters: Number of statistic_regions = ', &
1065                       statistic_regions+1
1066          PRINT*, '    Only 10 regions are allowed'
1067       ENDIF
1068       CALL local_stop
1069    ENDIF
1070    IF ( normalizing_region > statistic_regions  .OR. &
1071         normalizing_region < 0)  THEN
1072       IF ( myid == 0 )  THEN
1073          PRINT*, '+++ check_parameters: normalizing_region = ', &
1074                       normalizing_region, ' is unknown'
1075          PRINT*, '    Must be <= ', statistic_regions
1076       ENDIF
1077       CALL local_stop
1078    ENDIF
1079
1080!
1081!-- Set the default intervals for data output, if necessary
1082!-- NOTE: dt_dosp has already been set in package_parin
1083    IF ( dt_data_output /= 9999999.9 )  THEN
1084       IF ( dt_dopr           == 9999999.9 )  dt_dopr           = dt_data_output
1085       IF ( dt_dopts          == 9999999.9 )  dt_dopts          = dt_data_output
1086       IF ( dt_do2d_xy        == 9999999.9 )  dt_do2d_xy        = dt_data_output
1087       IF ( dt_do2d_xz        == 9999999.9 )  dt_do2d_xz        = dt_data_output
1088       IF ( dt_do2d_yz        == 9999999.9 )  dt_do2d_yz        = dt_data_output
1089       IF ( dt_do3d           == 9999999.9 )  dt_do3d           = dt_data_output
1090       IF ( dt_data_output_av == 9999999.9 )  dt_data_output_av = dt_data_output
1091    ENDIF
1092
1093!
1094!-- Set the default skip time intervals for data output, if necessary
1095    IF ( skip_time_dopr    == 9999999.9 ) &
1096                                       skip_time_dopr    = skip_time_data_output
1097    IF ( skip_time_dosp    == 9999999.9 ) &
1098                                       skip_time_dosp    = skip_time_data_output
1099    IF ( skip_time_do2d_xy == 9999999.9 ) &
1100                                       skip_time_do2d_xy = skip_time_data_output
1101    IF ( skip_time_do2d_xz == 9999999.9 ) &
1102                                       skip_time_do2d_xz = skip_time_data_output
1103    IF ( skip_time_do2d_yz == 9999999.9 ) &
1104                                       skip_time_do2d_yz = skip_time_data_output
1105    IF ( skip_time_do3d    == 9999999.9 ) &
1106                                       skip_time_do3d    = skip_time_data_output
1107    IF ( skip_time_data_output_av == 9999999.9 ) &
1108                                skip_time_data_output_av = skip_time_data_output
1109
1110!
1111!-- Check the average intervals (first for 3d-data, then for profiles and
1112!-- spectra)
1113    IF ( averaging_interval > dt_data_output_av )  THEN
1114       IF ( myid == 0 )  THEN
1115          PRINT*, '+++ check_parameters: average_interval=',              &
1116                       averaging_interval, ' must be <= dt_data_output=', &
1117                       dt_data_output
1118       ENDIF
1119       CALL local_stop
1120    ENDIF
1121
1122    IF ( averaging_interval_pr == 9999999.9 )  THEN
1123       averaging_interval_pr = averaging_interval
1124    ENDIF
1125
1126    IF ( averaging_interval_pr > dt_dopr )  THEN
1127       IF ( myid == 0 )  THEN
1128          PRINT*, '+++ check_parameters: averaging_interval_pr=', &
1129                       averaging_interval_pr, ' must be <= dt_dopr=', dt_dopr
1130       ENDIF
1131       CALL local_stop
1132    ENDIF
1133
1134    IF ( averaging_interval_sp == 9999999.9 )  THEN
1135       averaging_interval_sp = averaging_interval
1136    ENDIF
1137
1138    IF ( averaging_interval_sp > dt_dosp )  THEN
1139       IF ( myid == 0 )  THEN
1140          PRINT*, '+++ check_parameters: averaging_interval_sp=', &
1141                       averaging_interval_sp, ' must be <= dt_dosp=', dt_dosp
1142       ENDIF
1143       CALL local_stop
1144    ENDIF
1145
1146!
1147!-- Set the default interval for profiles entering the temporal average
1148    IF ( dt_averaging_input_pr == 9999999.9 )  THEN
1149       dt_averaging_input_pr = dt_averaging_input
1150    ENDIF
1151
1152!
1153!-- Set the default interval for the output of timeseries to a reasonable
1154!-- value (tries to minimize the number of calls of flow_statistics)
1155    IF ( dt_dots == 9999999.9 )  THEN
1156       IF ( averaging_interval_pr == 0.0 )  THEN
1157          dt_dots = MIN( dt_run_control, dt_dopr )
1158       ELSE
1159          dt_dots = MIN( dt_run_control, dt_averaging_input_pr )
1160       ENDIF
1161    ENDIF
1162
1163!
1164!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
1165    IF ( dt_averaging_input > averaging_interval )  THEN
1166       IF ( myid == 0 )  THEN
1167          PRINT*, '+++ check_parameters: dt_averaging_input=',                &
1168                       dt_averaging_input, ' must be <= averaging_interval=', &
1169                       averaging_interval
1170       ENDIF
1171       CALL local_stop
1172    ENDIF
1173
1174    IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
1175       IF ( myid == 0 )  THEN
1176          PRINT*, '+++ check_parameters: dt_averaging_input_pr=', &
1177                       dt_averaging_input_pr,                     &
1178                       ' must be <= averaging_interval_pr=',      &
1179                       averaging_interval_pr
1180       ENDIF
1181       CALL local_stop
1182    ENDIF
1183
1184!
1185!-- Set the default value for the integration interval of precipitation amount
1186    IF ( precipitation )  THEN
1187       IF ( precipitation_amount_interval == 9999999.9 )  THEN
1188          precipitation_amount_interval = dt_do2d_xy
1189       ELSE
1190          IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
1191             IF ( myid == 0 )  PRINT*, '+++ check_parameters: ',              &
1192                                       'precipitation_amount_interval =',     &
1193                                        precipitation_amount_interval,        &
1194                                       ' must not be larger than dt_do2d_xy', &
1195                                       ' = ', dt_do2d_xy   
1196       CALL local_stop
1197          ENDIF
1198       ENDIF
1199    ENDIF
1200
1201!
1202!-- Determine the number of output profiles and check whether they are
1203!-- permissible
1204    DO  WHILE ( data_output_pr(dopr_n+1) /= '          ' )
1205
1206       dopr_n = dopr_n + 1
1207       i = dopr_n
1208
1209!
1210!--    Determine internal profile number (for hom, homs)
1211!--    and store height levels
1212       SELECT CASE ( TRIM( data_output_pr(i) ) )
1213
1214          CASE ( 'u', '#u' )
1215             dopr_index(i) = 1
1216             hom(:,2,1,:)  = SPREAD( zu, 2, statistic_regions+1 )
1217             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1218                dopr_initial_index(i) = 5
1219                hom(:,2,5,:)          = SPREAD( zu, 2, statistic_regions+1 )
1220                data_output_pr(i)     = data_output_pr(i)(2:)
1221             ENDIF
1222
1223          CASE ( 'v', '#v' )
1224             dopr_index(i) = 2
1225             hom(:,2,2,:) = SPREAD( zu, 2, statistic_regions+1 )
1226             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1227                dopr_initial_index(i) = 6
1228                hom(:,2,6,:)          = SPREAD( zu, 2, statistic_regions+1 )
1229                data_output_pr(i)     = data_output_pr(i)(2:)
1230             ENDIF
1231
1232          CASE ( 'w' )
1233             dopr_index(i) = 3
1234             hom(:,2,3,:) = SPREAD( zw, 2, statistic_regions+1 )
1235
1236          CASE ( 'pt', '#pt' )
1237             IF ( .NOT. cloud_physics ) THEN
1238                dopr_index(i) = 4
1239                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
1240                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1241                   dopr_initial_index(i) = 7
1242                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
1243                   hom(nzb,2,7,:)        = 0.0    ! weil zu(nzb) negativ ist
1244                   data_output_pr(i)     = data_output_pr(i)(2:)
1245                ENDIF
1246             ELSE
1247                dopr_index(i) = 43
1248                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
1249                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1250                   dopr_initial_index(i) = 28
1251                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
1252                   hom(nzb,2,28,:)       = 0.0    ! weil zu(nzb) negativ ist
1253                   data_output_pr(i)     = data_output_pr(i)(2:)
1254                ENDIF
1255             ENDIF
1256
1257          CASE ( 'e' )
1258             dopr_index(i)  = 8
1259             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
1260             hom(nzb,2,8,:) = 0.0
1261
1262          CASE ( 'km', '#km' )
1263             dopr_index(i)  = 9
1264             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
1265             hom(nzb,2,9,:) = 0.0
1266             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1267                dopr_initial_index(i) = 23
1268                hom(:,2,23,:)         = hom(:,2,9,:)
1269                data_output_pr(i)     = data_output_pr(i)(2:)
1270             ENDIF
1271
1272          CASE ( 'kh', '#kh' )
1273             dopr_index(i)   = 10
1274             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
1275             hom(nzb,2,10,:) = 0.0
1276             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1277                dopr_initial_index(i) = 24
1278                hom(:,2,24,:)         = hom(:,2,10,:)
1279                data_output_pr(i)     = data_output_pr(i)(2:)
1280             ENDIF
1281
1282          CASE ( 'l', '#l' )
1283             dopr_index(i)   = 11
1284             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
1285             hom(nzb,2,11,:) = 0.0
1286             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1287                dopr_initial_index(i) = 25
1288                hom(:,2,25,:)         = hom(:,2,11,:)
1289                data_output_pr(i)     = data_output_pr(i)(2:)
1290             ENDIF
1291
1292          CASE ( 'w"u"' )
1293             dopr_index(i) = 12
1294             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
1295             IF ( prandtl_layer )  hom(nzb,2,12,:) = zu(1)
1296
1297          CASE ( 'w*u*' )
1298             dopr_index(i) = 13
1299             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
1300
1301          CASE ( 'w"v"' )
1302             dopr_index(i) = 14
1303             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
1304             IF ( prandtl_layer )  hom(nzb,2,14,:) = zu(1)
1305
1306          CASE ( 'w*v*' )
1307             dopr_index(i) = 15
1308             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
1309
1310          CASE ( 'w"pt"' )
1311             dopr_index(i) = 16
1312             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
1313
1314          CASE ( 'w*pt*' )
1315             dopr_index(i) = 17
1316             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
1317
1318          CASE ( 'wpt' )
1319             dopr_index(i) = 18
1320             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
1321
1322          CASE ( 'wu' )
1323             dopr_index(i) = 19
1324             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
1325             IF ( prandtl_layer )  hom(nzb,2,19,:) = zu(1)
1326
1327          CASE ( 'wv' )
1328             dopr_index(i) = 20
1329             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
1330             IF ( prandtl_layer )  hom(nzb,2,20,:) = zu(1)
1331
1332          CASE ( 'w*pt*BC' )
1333             dopr_index(i) = 21
1334             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
1335
1336          CASE ( 'wptBC' )
1337             dopr_index(i) = 22
1338             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
1339
1340          CASE ( 'u*2' )
1341             dopr_index(i) = 30
1342             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
1343
1344          CASE ( 'v*2' )
1345             dopr_index(i) = 31
1346             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
1347
1348          CASE ( 'w*2' )
1349             dopr_index(i) = 32
1350             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
1351
1352          CASE ( 'pt*2' )
1353             dopr_index(i) = 33
1354             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
1355
1356          CASE ( 'e*' )
1357             dopr_index(i) = 34
1358             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
1359
1360          CASE ( 'w*2pt*' )
1361             dopr_index(i) = 35
1362             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
1363
1364          CASE ( 'w*pt*2' )
1365             dopr_index(i) = 36
1366             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
1367
1368          CASE ( 'w*e*' )
1369             dopr_index(i) = 37
1370             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
1371
1372          CASE ( 'w*3' )
1373             dopr_index(i) = 38
1374             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
1375
1376          CASE ( 'Sw' )
1377             dopr_index(i) = 39
1378             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
1379
1380          CASE ( 'q', '#q' )
1381             IF ( .NOT. cloud_physics )  THEN
1382                IF ( myid == 0 )  THEN
1383                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1384                           data_output_pr(i),                          &
1385                           '    is not implemented for cloud_physics = FALSE'
1386                ENDIF
1387                CALL local_stop
1388             ELSE
1389                dopr_index(i) = 41
1390                hom(:,2,41,:)  = SPREAD( zu, 2, statistic_regions+1 )
1391                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1392                   dopr_initial_index(i) = 26
1393                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
1394                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
1395                   data_output_pr(i)     = data_output_pr(i)(2:)
1396                ENDIF
1397             ENDIF
1398
1399          CASE ( 's', '#s' )
1400             IF ( .NOT. passive_scalar )  THEN
1401                IF ( myid == 0 )  THEN
1402                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1403                           data_output_pr(i),                          &
1404                           '    is not implemented for passive_scalar = FALSE'
1405                ENDIF
1406                CALL local_stop
1407             ELSE
1408                dopr_index(i) = 41
1409                hom(:,2,41,:)  = SPREAD( zu, 2, statistic_regions+1 )
1410                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1411                   dopr_initial_index(i) = 26
1412                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
1413                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
1414                   data_output_pr(i)     = data_output_pr(i)(2:)
1415                ENDIF
1416             ENDIF
1417
1418          CASE ( 'qv', '#qv' )
1419             IF ( .NOT. cloud_physics ) THEN
1420                dopr_index(i) = 41
1421                hom(:,2,41,:)  = SPREAD( zu, 2, statistic_regions+1 )
1422                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1423                   dopr_initial_index(i) = 26
1424                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
1425                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
1426                   data_output_pr(i)     = data_output_pr(i)(2:)
1427                ENDIF
1428             ELSE
1429                dopr_index(i) = 42
1430                hom(:,2,42,:)  = SPREAD( zu, 2, statistic_regions+1 )
1431                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1432                   dopr_initial_index(i) = 27
1433                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
1434                   hom(nzb,2,27,:)       = 0.0    ! weil zu(nzb) negativ ist
1435                   data_output_pr(i)     = data_output_pr(i)(2:)
1436                ENDIF
1437             ENDIF
1438
1439          CASE ( 'lpt', '#lpt' )
1440             IF ( .NOT. cloud_physics ) THEN
1441                IF ( myid == 0 )  THEN
1442                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1443                           data_output_pr(i),                          &
1444                           '    is not implemented for cloud_physics = FALSE'
1445                ENDIF
1446                CALL local_stop
1447             ELSE
1448                dopr_index(i) = 4
1449                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
1450                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1451                   dopr_initial_index(i) = 7
1452                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
1453                   hom(nzb,2,7,:)        = 0.0    ! weil zu(nzb) negativ ist
1454                   data_output_pr(i)     = data_output_pr(i)(2:)
1455                ENDIF
1456             ENDIF
1457
1458          CASE ( 'vpt', '#vpt' )
1459             dopr_index(i) = 44
1460             hom(:,2,44,:)  = SPREAD( zu, 2, statistic_regions+1 )
1461             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1462                dopr_initial_index(i) = 29
1463                hom(:,2,29,:)         = SPREAD( zu, 2, statistic_regions+1 )
1464                hom(nzb,2,29,:)       = 0.0    ! weil zu(nzb) negativ ist
1465                data_output_pr(i)     = data_output_pr(i)(2:)
1466             ENDIF
1467
1468          CASE ( 'w"vpt"' )
1469             dopr_index(i) = 45
1470             hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
1471
1472          CASE ( 'w*vpt*' )
1473             dopr_index(i) = 46
1474             hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
1475
1476          CASE ( 'wvpt' )
1477             dopr_index(i) = 47
1478             hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
1479
1480          CASE ( 'w"q"' )
1481             IF ( .NOT. cloud_physics ) THEN
1482                IF ( myid == 0 )  THEN
1483                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1484                           data_output_pr(i),                          &
1485                           '    is not implemented for cloud_physics = FALSE'
1486                ENDIF
1487                CALL local_stop
1488             ELSE
1489                dopr_index(i) = 48
1490                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
1491             ENDIF
1492
1493          CASE ( 'w*q*' )
1494             IF ( .NOT. cloud_physics ) THEN
1495                IF ( myid == 0 )  THEN
1496                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1497                           data_output_pr(i),                          &
1498                           '    is not implemented for cloud_physics = FALSE'
1499                ENDIF
1500                CALL local_stop
1501             ELSE
1502                dopr_index(i) = 49
1503                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
1504             ENDIF
1505
1506          CASE ( 'wq' )
1507             IF ( .NOT. cloud_physics ) THEN
1508                IF ( myid == 0 )  THEN
1509                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1510                           data_output_pr(i),                          &
1511                           '    is not implemented for cloud_physics = FALSE'
1512                ENDIF
1513                CALL local_stop
1514             ELSE
1515                dopr_index(i) = 50
1516                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
1517             ENDIF
1518
1519          CASE ( 'w"s"' )
1520             IF ( .NOT. passive_scalar ) THEN
1521                IF ( myid == 0 )  THEN
1522                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1523                           data_output_pr(i),                          &
1524                           '    is not implemented for passive_scalar = FALSE'
1525                ENDIF
1526                CALL local_stop
1527             ELSE
1528                dopr_index(i) = 48
1529                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
1530             ENDIF
1531
1532          CASE ( 'w*s*' )
1533             IF ( .NOT. passive_scalar ) THEN
1534                IF ( myid == 0 )  THEN
1535                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1536                           data_output_pr(i),                          &
1537                           '    is not implemented for passive_scalar = FALSE'
1538                ENDIF
1539                CALL local_stop
1540             ELSE
1541                dopr_index(i) = 49
1542                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
1543             ENDIF
1544
1545          CASE ( 'ws' )
1546             IF ( .NOT. passive_scalar ) THEN
1547                IF ( myid == 0 )  THEN
1548                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1549                           data_output_pr(i),                          &
1550                           '    is not implemented for passive_scalar = FALSE'
1551                ENDIF
1552                CALL local_stop
1553             ELSE
1554                dopr_index(i) = 50
1555                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
1556             ENDIF
1557
1558          CASE ( 'w"qv"' )
1559             IF ( humidity  .AND.  .NOT. cloud_physics ) &
1560             THEN
1561                dopr_index(i) = 48
1562                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
1563             ELSEIF( humidity .AND. cloud_physics ) THEN
1564                dopr_index(i) = 51
1565                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
1566             ELSE
1567                IF ( myid == 0 )  THEN
1568                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1569                           data_output_pr(i),                          &
1570                           '    is not implemented for cloud_physics = FALSE', &
1571                           '    and                    humidity      = FALSE'
1572                ENDIF
1573                CALL local_stop                   
1574             ENDIF
1575
1576          CASE ( 'w*qv*' )
1577             IF ( humidity  .AND.  .NOT. cloud_physics ) &
1578             THEN
1579                dopr_index(i) = 49
1580                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
1581             ELSEIF( humidity .AND. cloud_physics ) THEN
1582                dopr_index(i) = 52
1583                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
1584             ELSE
1585                IF ( myid == 0 )  THEN
1586                   PRINT*, '+++ check_parameters:  data_output_pr = ',         &
1587                           data_output_pr(i),                                  &
1588                           '    is not implemented for cloud_physics = FALSE', &
1589                           '                       and humidity      = FALSE'
1590                ENDIF
1591                CALL local_stop                   
1592             ENDIF
1593
1594          CASE ( 'wqv' )
1595             IF ( humidity  .AND.  .NOT. cloud_physics ) &
1596             THEN
1597                dopr_index(i) = 50
1598                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
1599             ELSEIF( humidity .AND. cloud_physics ) THEN
1600                dopr_index(i) = 53
1601                hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
1602             ELSE
1603                IF ( myid == 0 )  THEN
1604                   PRINT*, '+++ check_parameters:  data_output_pr = ',         &
1605                           data_output_pr(i),                                  &
1606                           '    is not implemented for cloud_physics = FALSE', &
1607                           '                       and humidity      = FALSE'
1608                ENDIF
1609                CALL local_stop                   
1610             ENDIF
1611
1612          CASE ( 'ql' )
1613             IF ( .NOT. cloud_physics  .AND.  .NOT. cloud_droplets )  THEN
1614                IF ( myid == 0 )  THEN
1615                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1616                           data_output_pr(i),                          &
1617                           '    is not implemented for cloud_physics = FALSE'
1618                ENDIF
1619                CALL local_stop
1620             ELSE
1621                dopr_index(i) = 54
1622                hom(:,2,54,:)  = SPREAD( zu, 2, statistic_regions+1 )
1623             ENDIF
1624
1625          CASE ( 'w*u*u*/dz' )
1626             dopr_index(i) = 55
1627             hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
1628
1629          CASE ( 'w*p*/dz' )
1630             dopr_index(i) = 56
1631             hom(:,2,56,:) = SPREAD( zu, 2, statistic_regions+1 )
1632
1633          CASE ( 'w"e/dz' )
1634             dopr_index(i) = 57
1635             hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
1636
1637          CASE ( 'u"pt"' )
1638             dopr_index(i) = 58
1639             hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
1640
1641          CASE ( 'u*pt*' )
1642             dopr_index(i) = 59
1643             hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
1644
1645          CASE ( 'upt_t' )
1646             dopr_index(i) = 60
1647             hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
1648
1649          CASE ( 'v"pt"' )
1650             dopr_index(i) = 61
1651             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
1652             
1653          CASE ( 'v*pt*' )
1654             dopr_index(i) = 62
1655             hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
1656
1657          CASE ( 'vpt_t' )
1658             dopr_index(i) = 63
1659             hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
1660
1661
1662          CASE DEFAULT
1663             IF ( myid == 0 )  THEN
1664                PRINT*, '+++ check_parameters:  unknown output profile:  ', &
1665                        'data_output_pr = ', data_output_pr(i)
1666             ENDIF
1667             CALL local_stop
1668
1669       END SELECT
1670!
1671!--    Check to which of the predefined coordinate systems the profile belongs
1672       DO  k = 1, crmax
1673          IF ( INDEX( cross_profiles(k), ' '//TRIM( data_output_pr(i) )//' ' ) &
1674               /=0 ) &
1675          THEN
1676             dopr_crossindex(i) = k
1677             EXIT
1678          ENDIF
1679       ENDDO
1680!
1681!--    Generate the text for the labels of the PROFIL output file. "-characters
1682!--    must be substituted, otherwise PROFIL would interpret them as TeX
1683!--    control characters
1684       dopr_label(i) = data_output_pr(i)
1685       position = INDEX( dopr_label(i) , '"' )
1686       DO WHILE ( position /= 0 )
1687          dopr_label(i)(position:position) = ''''
1688          position = INDEX( dopr_label(i) , '"' )
1689       ENDDO
1690
1691    ENDDO
1692
1693!
1694!-- y-value range of the coordinate system (PROFIL).
1695!-- x-value range determined in plot_1d.
1696    cross_uymin = 0.0
1697    IF ( z_max_do1d == -1.0 )  THEN
1698       cross_uymax = zu(nzt+1)
1699    ELSEIF ( z_max_do1d < zu(nzb+1)  .OR.  z_max_do1d > zu(nzt+1) )  THEN
1700       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  z_max_do1d=',  &
1701                                 z_max_do1d,' must be >= ', zu(nzb+1),  &
1702                                 ' or <= ', zu(nzt+1)
1703       CALL local_stop
1704    ELSE
1705       cross_uymax = z_max_do1d
1706    ENDIF
1707
1708!
1709!-- Check whether the chosen normalizing factor for the coordinate systems is
1710!-- permissible
1711    DO  i = 1, crmax
1712       SELECT CASE ( TRIM( cross_normalized_x(i) ) )  ! TRIM required on IBM
1713
1714          CASE ( '', 'wpt0', 'ws2', 'tsw2', 'ws3', 'ws2tsw', 'wstsw2' )
1715             j = 0
1716
1717          CASE DEFAULT
1718             IF ( myid == 0 )  THEN
1719                PRINT*, '+++ check_parameters: unknown normalize method'
1720                PRINT*, '    cross_normalized_x="',cross_normalized_x(i),'"'
1721             ENDIF
1722             CALL local_stop
1723
1724       END SELECT
1725       SELECT CASE ( TRIM( cross_normalized_y(i) ) )  ! TRIM required on IBM
1726
1727          CASE ( '', 'z_i' )
1728             j = 0
1729
1730          CASE DEFAULT
1731             IF ( myid == 0 )  THEN
1732                PRINT*, '+++ check_parameters: unknown normalize method'
1733                PRINT*, '    cross_normalized_y="',cross_normalized_y(i),'"'
1734             ENDIF
1735             CALL local_stop
1736
1737       END SELECT
1738    ENDDO
1739!
1740!-- Check normalized y-value range of the coordinate system (PROFIL)
1741    IF ( z_max_do1d_normalized /= -1.0  .AND.  z_max_do1d_normalized <= 0.0 ) &
1742    THEN
1743       IF ( myid == 0 )  PRINT*,'+++ check_parameters:  z_max_do1d_normalize', &
1744                                'd=', z_max_do1d_normalized, ' must be >= 0.0'
1745       CALL local_stop
1746    ENDIF
1747
1748
1749!
1750!-- Append user-defined data output variables to the standard data output
1751    IF ( data_output_user(1) /= ' ' )  THEN
1752       i = 1
1753       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
1754          i = i + 1
1755       ENDDO
1756       j = 1
1757       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 100 )
1758          IF ( i > 100 )  THEN
1759             IF ( myid == 0 )  THEN
1760                PRINT*, '+++ check_parameters: number of output quantitities', &
1761                             ' given by data_output and data_output_user'
1762                PRINT*, '                      exceeds the limit of 100'
1763             ENDIF
1764             CALL local_stop
1765          ENDIF
1766          data_output(i) = data_output_user(j)
1767          i = i + 1
1768          j = j + 1
1769       ENDDO
1770    ENDIF
1771
1772!
1773!-- Check and set steering parameters for 2d/3d data output and averaging
1774    i   = 1
1775    DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
1776!
1777!--    Check for data averaging
1778       ilen = LEN_TRIM( data_output(i) )
1779       j = 0                                                 ! no data averaging
1780       IF ( ilen > 3 )  THEN
1781          IF ( data_output(i)(ilen-2:ilen) == '_av' )  THEN
1782             j = 1                                           ! data averaging
1783             data_output(i) = data_output(i)(1:ilen-3)
1784          ENDIF
1785       ENDIF
1786!
1787!--    Check for cross section or volume data
1788       ilen = LEN_TRIM( data_output(i) )
1789       k = 0                                                   ! 3d data
1790       var = data_output(i)(1:ilen)
1791       IF ( ilen > 3 )  THEN
1792          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR. &
1793               data_output(i)(ilen-2:ilen) == '_xz'  .OR. &
1794               data_output(i)(ilen-2:ilen) == '_yz' )  THEN
1795             k = 1                                             ! 2d data
1796             var = data_output(i)(1:ilen-3)
1797          ENDIF
1798       ENDIF
1799!
1800!--    Check for allowed value and set units
1801       SELECT CASE ( TRIM( var ) )
1802
1803          CASE ( 'e' )
1804             IF ( constant_diffusion )  THEN
1805                IF ( myid == 0 )  THEN
1806                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1807                                '" requires constant_diffusion = .FALSE.'
1808                ENDIF
1809                CALL local_stop
1810             ENDIF
1811             unit = 'm2/s2'
1812
1813          CASE ( 'pc', 'pr' )
1814             IF ( .NOT. particle_advection )  THEN
1815                IF ( myid == 0 )  THEN
1816                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1817                                '" requires particle package'
1818                   PRINT*, '                      (mrun-option "-p particles")'
1819                ENDIF
1820                CALL local_stop
1821             ENDIF
1822             IF ( TRIM( var ) == 'pc' )  unit = 'number'
1823             IF ( TRIM( var ) == 'pr' )  unit = 'm'
1824
1825          CASE ( 'q', 'vpt' )
1826             IF ( .NOT. humidity )  THEN
1827                IF ( myid == 0 )  THEN
1828                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1829                                '" requires humidity = .TRUE.'
1830                ENDIF
1831                CALL local_stop
1832             ENDIF
1833             IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
1834             IF ( TRIM( var ) == 'vpt' )  unit = 'K'
1835
1836          CASE ( 'ql' )
1837             IF ( .NOT. ( cloud_physics  .OR.  cloud_droplets ) )  THEN
1838                IF ( myid == 0 )  THEN
1839                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1840                                '" requires cloud_physics = .TRUE.'
1841                   PRINT*, '                      or cloud_droplets = .TRUE.'
1842                ENDIF
1843                CALL local_stop
1844             ENDIF
1845             unit = 'kg/kg'
1846
1847          CASE ( 'ql_c', 'ql_v', 'ql_vp' )
1848             IF ( .NOT. cloud_droplets )  THEN
1849                IF ( myid == 0 )  THEN
1850                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1851                                '" requires cloud_droplets = .TRUE.'
1852                ENDIF
1853                CALL local_stop
1854             ENDIF
1855             IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
1856             IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
1857             IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
1858
1859          CASE ( 'qv' )
1860             IF ( .NOT. cloud_physics )  THEN
1861                IF ( myid == 0 )  THEN
1862                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1863                                '" requires cloud_physics = .TRUE.'
1864                ENDIF
1865                CALL local_stop
1866             ENDIF
1867             unit = 'kg/kg'
1868
1869          CASE ( 's' )
1870             IF ( .NOT. passive_scalar )  THEN
1871                IF ( myid == 0 )  THEN
1872                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1873                                '" requires passive_scalar = .TRUE.'
1874                ENDIF
1875                CALL local_stop
1876             ENDIF
1877             unit = 'conc'
1878
1879          CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*' )
1880             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
1881                IF ( myid == 0 )  THEN
1882                   PRINT*, '+++ check_parameters:  illegal value for data_',&
1883                                'output: "', TRIM( var ), '" is only allowed'
1884                   PRINT*, '                       for horizontal cross section'
1885                ENDIF
1886                CALL local_stop
1887             ENDIF
1888             IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
1889                IF ( myid == 0 )  THEN
1890                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1891                                '" requires cloud_physics = .TRUE.'
1892                ENDIF
1893                CALL local_stop
1894             ENDIF
1895             IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
1896                IF ( myid == 0 )  THEN
1897                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1898                                '" requires precipitation = .TRUE.'
1899                ENDIF
1900                CALL local_stop
1901             ENDIF
1902             IF ( TRIM( var ) == 'pra*'  .AND.  j == 1 )  THEN
1903                IF ( myid == 0 )  THEN
1904                   PRINT*, '+++ check_parameters: temporal averaging of ', &
1905                           ' precipitation amount "', TRIM( var ),         &
1906                           '" not possible' 
1907                ENDIF
1908                CALL local_stop
1909             ENDIF
1910             IF ( TRIM( var ) == 'prr*'  .AND.  .NOT. precipitation )  THEN
1911                IF ( myid == 0 )  THEN
1912                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1913                                '" requires precipitation = .TRUE.'
1914                ENDIF
1915                CALL local_stop
1916             ENDIF
1917
1918
1919             IF ( TRIM( var ) == 'u*'   )  unit = 'm/s'
1920             IF ( TRIM( var ) == 't*'   )  unit = 'K'
1921             IF ( TRIM( var ) == 'lwp*' )  unit = 'kg/kg*m'
1922             IF ( TRIM( var ) == 'pra*' )  unit = 'mm'
1923             IF ( TRIM( var ) == 'prr*' )  unit = 'mm/s'
1924             IF ( TRIM( var ) == 'z0*'  )  unit = 'm'
1925
1926          CASE ( 'p', 'pt', 'u', 'v', 'w' )
1927             IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
1928             IF ( TRIM( var ) == 'pt' )  unit = 'K'
1929             IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
1930             IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
1931             IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
1932             CONTINUE
1933
1934          CASE DEFAULT
1935             CALL user_check_data_output( var, unit )
1936
1937             IF ( unit == 'illegal' )  THEN
1938                IF ( myid == 0 )  THEN
1939                   IF ( data_output_user(1) /= ' ' )  THEN
1940                      PRINT*, '+++ check_parameters:  illegal value for data_',&
1941                                   'output or data_output_user: "',            &
1942                                   TRIM( data_output(i) ), '"'
1943                   ELSE
1944                      PRINT*, '+++ check_parameters:  illegal value for data_',&
1945                                   'output: "', TRIM( data_output(i) ), '"'
1946                   ENDIF
1947                ENDIF
1948                CALL local_stop
1949             ENDIF
1950
1951       END SELECT
1952!
1953!--    Set the internal steering parameters appropriately
1954       IF ( k == 0 )  THEN
1955          do3d_no(j)              = do3d_no(j) + 1
1956          do3d(j,do3d_no(j))      = data_output(i)
1957          do3d_unit(j,do3d_no(j)) = unit
1958       ELSE
1959          do2d_no(j)              = do2d_no(j) + 1
1960          do2d(j,do2d_no(j))      = data_output(i)
1961          do2d_unit(j,do2d_no(j)) = unit
1962          IF ( data_output(i)(ilen-2:ilen) == '_xy' )  THEN
1963             data_output_xy(j) = .TRUE.
1964          ENDIF
1965          IF ( data_output(i)(ilen-2:ilen) == '_xz' )  THEN
1966             data_output_xz(j) = .TRUE.
1967          ENDIF
1968          IF ( data_output(i)(ilen-2:ilen) == '_yz' )  THEN
1969             data_output_yz(j) = .TRUE.
1970          ENDIF
1971       ENDIF
1972
1973       IF ( j == 1 )  THEN
1974!
1975!--       Check, if variable is already subject to averaging
1976          found = .FALSE.
1977          DO  k = 1, doav_n
1978             IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
1979          ENDDO
1980
1981          IF ( .NOT. found )  THEN
1982             doav_n = doav_n + 1
1983             doav(doav_n) = var
1984          ENDIF
1985       ENDIF
1986
1987       i = i + 1
1988    ENDDO
1989
1990!
1991!-- Store sectional planes in one shared array
1992    section(:,1) = section_xy
1993    section(:,2) = section_xz
1994    section(:,3) = section_yz
1995
1996!
1997!-- Upper plot limit (grid point value) for 1D profiles
1998    IF ( z_max_do1d == -1.0 )  THEN
1999       nz_do1d = nzt+1
2000    ELSE
2001       DO  k = nzb+1, nzt+1
2002          nz_do1d = k
2003          IF ( zw(k) > z_max_do1d )  EXIT
2004       ENDDO
2005    ENDIF
2006
2007!
2008!-- Upper plot limit for 2D vertical sections
2009    IF ( z_max_do2d == -1.0 )  z_max_do2d = zu(nzt)
2010    IF ( z_max_do2d < zu(nzb+1)  .OR.  z_max_do2d > zu(nzt) )  THEN
2011       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  z_max_do2d=',        &
2012                                 z_max_do2d, ' must be >= ', zu(nzb+1),       &
2013                                 '(zu(nzb+1)) and <= ', zu(nzt), ' (zu(nzt))'
2014       CALL local_stop
2015    ENDIF
2016
2017!
2018!-- Upper plot limit for 3D arrays
2019    IF ( nz_do3d == -9999 )  nz_do3d = nzt + 1
2020
2021!
2022!-- Determine and check accuracy for compressed 3D plot output
2023    IF ( do3d_compress )  THEN
2024!
2025!--    Compression only permissible on T3E machines
2026       IF ( host(1:3) /= 't3e' )  THEN
2027          IF ( myid == 0 )  THEN
2028             PRINT*, '+++ check_parameters: do3d_compress = .TRUE. not allow', &
2029                          'ed on host "', TRIM( host ), '"'
2030          ENDIF
2031          CALL local_stop
2032       ENDIF
2033
2034       i = 1
2035       DO  WHILE ( do3d_comp_prec(i) /= ' ' )
2036
2037          ilen = LEN_TRIM( do3d_comp_prec(i) )
2038          IF ( LLT( do3d_comp_prec(i)(ilen:ilen), '0' ) .OR. &
2039               LGT( do3d_comp_prec(i)(ilen:ilen), '9' ) )  THEN
2040             IF ( myid == 0 )  THEN
2041                PRINT*, '+++ check_parameters: illegal precision: ', &
2042                        'do3d_comp_prec(', i, ')="', TRIM(do3d_comp_prec(i)),'"'
2043             ENDIF
2044             CALL local_stop
2045          ENDIF
2046
2047          prec = IACHAR( do3d_comp_prec(i)(ilen:ilen) ) - IACHAR( '0' )
2048          var = do3d_comp_prec(i)(1:ilen-1)
2049
2050          SELECT CASE ( var )
2051
2052             CASE ( 'u' )
2053                j = 1
2054             CASE ( 'v' )
2055                j = 2
2056             CASE ( 'w' )
2057                j = 3
2058             CASE ( 'p' )
2059                j = 4
2060             CASE ( 'pt' )
2061                j = 5
2062
2063             CASE DEFAULT
2064                IF ( myid == 0 )  THEN
2065                   PRINT*, '+++ check_parameters: unknown variable in ', &
2066                               'assignment'
2067                   PRINT*, '    do3d_comp_prec(', i, ')="', &
2068                           TRIM( do3d_comp_prec(i) ),'"'
2069                ENDIF
2070                CALL local_stop               
2071
2072          END SELECT
2073
2074          plot_3d_precision(j)%precision = prec
2075          i = i + 1
2076
2077       ENDDO
2078    ENDIF
2079
2080!
2081!-- Check the data output format(s)
2082    IF ( data_output_format(1) == ' ' )  THEN
2083!
2084!--    Default value
2085       netcdf_output = .TRUE.
2086    ELSE
2087       i = 1
2088       DO  WHILE ( data_output_format(i) /= ' ' )
2089
2090          SELECT CASE ( data_output_format(i) )
2091
2092             CASE ( 'netcdf' )
2093                netcdf_output = .TRUE.
2094             CASE ( 'iso2d' )
2095                iso2d_output  = .TRUE.
2096             CASE ( 'profil' )
2097                profil_output = .TRUE.
2098             CASE ( 'avs' )
2099                avs_output    = .TRUE.
2100
2101             CASE DEFAULT
2102                IF ( myid == 0 )  THEN
2103                   PRINT*, '+++ check_parameters:'
2104                   PRINT*, '    unknown value for data_output_format "', &
2105                                TRIM( data_output_format(i) ),'"'
2106                ENDIF
2107                CALL local_stop               
2108
2109          END SELECT
2110
2111          i = i + 1
2112          IF ( i > 10 )  EXIT
2113
2114       ENDDO
2115
2116    ENDIF
2117
2118!
2119!-- Check netcdf precison
2120    ldum = .FALSE.
2121    CALL define_netcdf_header( 'ch', ldum, 0 )
2122
2123!
2124!-- Check, whether a constant diffusion coefficient shall be used
2125    IF ( km_constant /= -1.0 )  THEN
2126       IF ( km_constant < 0.0 )  THEN
2127          IF ( myid == 0 )  PRINT*, '+++ check_parameters:  km_constant=', &
2128                                    km_constant, ' < 0.0'
2129          CALL local_stop
2130       ELSE
2131          IF ( prandtl_number < 0.0 )  THEN
2132             IF ( myid == 0 )  PRINT*,'+++ check_parameters:  prandtl_number=',&
2133                                      prandtl_number, ' < 0.0'
2134             CALL local_stop
2135          ENDIF
2136          constant_diffusion = .TRUE.
2137
2138          IF ( prandtl_layer )  THEN
2139             IF ( myid == 0 )  PRINT*, '+++ check_parameters:  prandtl_layer ',&
2140                                       'is not allowed with fixed value of km'
2141             CALL local_stop
2142          ENDIF
2143       ENDIF
2144    ENDIF
2145
2146!
2147!-- In case of non-cyclic lateral boundaries, set the default maximum value
2148!-- for the horizontal diffusivity used within the outflow damping layer,
2149!-- and check/set the width of the damping layer
2150    IF ( bc_lr /= 'cyclic' ) THEN
2151       IF ( km_damp_max == -1.0 )  THEN
2152          km_damp_max = 0.5 * dx
2153       ENDIF
2154       IF ( outflow_damping_width == -1.0 )  THEN
2155          outflow_damping_width = MIN( 20, nx/2 )
2156       ENDIF
2157       IF ( outflow_damping_width <= 0  .OR.  outflow_damping_width > nx )  THEN
2158          IF ( myid == 0 )  PRINT*, '+++ check_parameters: outflow_damping w',&
2159                                    'idth out of range'
2160          CALL local_stop
2161       ENDIF
2162    ENDIF
2163
2164    IF ( bc_ns /= 'cyclic' )  THEN
2165       IF ( km_damp_max == -1.0 )  THEN
2166          km_damp_max = 0.5 * dy
2167       ENDIF
2168       IF ( outflow_damping_width == -1.0 )  THEN
2169          outflow_damping_width = MIN( 20, ny/2 )
2170       ENDIF
2171       IF ( outflow_damping_width <= 0  .OR.  outflow_damping_width > ny )  THEN
2172          IF ( myid == 0 )  PRINT*, '+++ check_parameters: outflow_damping w',&
2173                                    'idth out of range'
2174          CALL local_stop
2175       ENDIF
2176    ENDIF
2177
2178!
2179!-- Check value range for rif
2180    IF ( rif_min >= rif_max )  THEN
2181       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  rif_min=', rif_min, &
2182                                 ' must be less than rif_max=', rif_max
2183       CALL local_stop
2184    ENDIF
2185
2186!
2187!-- Determine upper and lower hight level indices for random perturbations
2188    IF ( disturbance_level_b == -1.0 )  THEN
2189       disturbance_level_b     = zu(nzb+3)
2190       disturbance_level_ind_b = nzb + 3
2191    ELSEIF ( disturbance_level_b < zu(3) )  THEN
2192       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_b=',&
2193                                 disturbance_level_b, ' must be >= ',zu(3),    &
2194                                 '(zu(3))'
2195       CALL local_stop
2196    ELSEIF ( disturbance_level_b > zu(nzt-2) )  THEN
2197       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_b=',&
2198                                 disturbance_level_b, ' must be <= ',zu(nzt-2),&
2199                                 '(zu(nzt-2))'
2200       CALL local_stop
2201    ELSE
2202       DO  k = 3, nzt-2
2203          IF ( disturbance_level_b <= zu(k) )  THEN
2204             disturbance_level_ind_b = k
2205             EXIT
2206          ENDIF
2207       ENDDO
2208    ENDIF
2209
2210    IF ( disturbance_level_t == -1.0 )  THEN
2211       disturbance_level_t     = zu(nzt/3)
2212       disturbance_level_ind_t = nzt / 3
2213    ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
2214       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_t=',&
2215                                 disturbance_level_t, ' must be <= ',zu(nzt-2),&
2216                                 '(zu(nzt-2))'
2217       CALL local_stop
2218    ELSEIF ( disturbance_level_t < disturbance_level_b )  THEN
2219       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_t=',&
2220                                 disturbance_level_t, ' must be >= ',          &
2221                                 'disturbance_level_b=', disturbance_level_b
2222       CALL local_stop
2223    ELSE
2224       DO  k = 3, nzt-2
2225          IF ( disturbance_level_t <= zu(k) )  THEN
2226             disturbance_level_ind_t = k
2227             EXIT
2228          ENDIF
2229       ENDDO
2230    ENDIF
2231
2232!
2233!-- Check again whether the levels determined this way are ok.
2234!-- Error may occur at automatic determination and too few grid points in
2235!-- z-direction.
2236    IF ( disturbance_level_ind_t < disturbance_level_ind_b )  THEN
2237       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  ',                &
2238                                 'disturbance_level_ind_t=',               &
2239                                 disturbance_level_ind_t, ' must be >= ',  &
2240                                 'disturbance_level_ind_b=',               &
2241                                 disturbance_level_ind_b
2242       CALL local_stop
2243    ENDIF
2244
2245!
2246!-- Determine the horizontal index range for random perturbations.
2247!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed
2248!-- near the inflow and the perturbation area is further limited to ...(1)
2249!-- after the initial phase of the flow.
2250    dist_nxl = 0;  dist_nxr = nx
2251    dist_nys = 0;  dist_nyn = ny
2252    IF ( bc_lr /= 'cyclic' )  THEN
2253       IF ( inflow_disturbance_begin == -1 )  THEN
2254          inflow_disturbance_begin = MIN( 10, nx/2 )
2255       ENDIF
2256       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )&
2257       THEN
2258          IF ( myid == 0 )  PRINT*, '+++ check_parameters: inflow_disturbance',&
2259                                    '_begin out of range'
2260          CALL local_stop
2261       ENDIF
2262       IF ( inflow_disturbance_end == -1 )  THEN
2263          inflow_disturbance_end = MIN( 100, 3*nx/4 )
2264       ENDIF
2265       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )    &
2266       THEN
2267          IF ( myid == 0 )  PRINT*, '+++ check_parameters: inflow_disturbance',&
2268                                    '_end out of range'
2269          CALL local_stop
2270       ENDIF
2271    ELSEIF ( bc_ns /= 'cyclic' )  THEN
2272       IF ( inflow_disturbance_begin == -1 )  THEN
2273          inflow_disturbance_begin = MIN( 10, ny/2 )
2274       ENDIF
2275       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )&
2276       THEN
2277          IF ( myid == 0 )  PRINT*, '+++ check_parameters: inflow_disturbance',&
2278                                    '_begin out of range'
2279          CALL local_stop
2280       ENDIF
2281       IF ( inflow_disturbance_end == -1 )  THEN
2282          inflow_disturbance_end = MIN( 100, 3*ny/4 )
2283       ENDIF
2284       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )    &
2285       THEN
2286          IF ( myid == 0 )  PRINT*, '+++ check_parameters: inflow_disturbance',&
2287                                    '_end out of range'
2288          CALL local_stop
2289       ENDIF
2290    ENDIF
2291
2292    IF ( bc_lr == 'radiation/dirichlet' )  THEN
2293       dist_nxr    = nx - inflow_disturbance_begin
2294       dist_nxl(1) = nx - inflow_disturbance_end
2295    ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
2296       dist_nxl    = inflow_disturbance_begin
2297       dist_nxr(1) = inflow_disturbance_end
2298    ENDIF
2299    IF ( bc_ns == 'dirichlet/radiation' )  THEN
2300       dist_nyn    = ny - inflow_disturbance_begin
2301       dist_nys(1) = ny - inflow_disturbance_end
2302    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
2303       dist_nys    = inflow_disturbance_begin
2304       dist_nyn(1) = inflow_disturbance_end
2305    ENDIF
2306
2307!
2308!-- Check random generator
2309    IF ( random_generator /= 'system-specific'  .AND. &
2310         random_generator /= 'numerical-recipes' )  THEN
2311       IF ( myid == 0 )  THEN
2312          PRINT*, '+++ check_parameters:'
2313          PRINT*, '    unknown random generator: random_generator=', &
2314                  random_generator
2315       ENDIF
2316       CALL local_stop
2317    ENDIF
2318
2319!
2320!-- Determine damping level index for 1D model
2321    IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
2322       IF ( damp_level_1d == -1.0 )  THEN
2323          damp_level_1d     = zu(nzt+1)
2324          damp_level_ind_1d = nzt + 1
2325       ELSEIF ( damp_level_1d < 0.0  .OR.  damp_level_1d > zu(nzt+1) )  THEN
2326          IF ( myid == 0 )  PRINT*, '+++ check_parameters:  damp_level_1d=', &
2327                                    damp_level_1d, ' must be > 0.0 and < ',  &
2328                                    'zu(nzt+1)', zu(nzt+1)
2329          CALL local_stop
2330       ELSE
2331          DO  k = 1, nzt+1
2332             IF ( damp_level_1d <= zu(k) )  THEN
2333                damp_level_ind_1d = k
2334                EXIT
2335             ENDIF
2336          ENDDO
2337       ENDIF
2338    ENDIF
2339!
2340!-- Check some other 1d-model parameters
2341    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND. &
2342         TRIM( mixing_length_1d ) /= 'blackadar' )  THEN
2343       IF ( myid == 0 )  PRINT*, '+++ check_parameters: mixing_length_1d = "', &
2344                                 TRIM( mixing_length_1d ), '" is unknown'
2345       CALL local_stop
2346    ENDIF
2347    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND. &
2348         TRIM( dissipation_1d ) /= 'detering' )  THEN
2349       IF ( myid == 0 )  PRINT*, '+++ check_parameters: dissipation_1d = "', &
2350                                 TRIM( dissipation_1d ), '" is unknown'
2351       CALL local_stop
2352    ENDIF
2353
2354!
2355!-- Set time for the next user defined restart (time_restart is the
2356!-- internal parameter for steering restart events)
2357    IF ( restart_time /= 9999999.9 )  THEN
2358       IF ( restart_time > simulated_time )  time_restart = restart_time
2359    ELSE
2360!
2361!--    In case of a restart run, set internal parameter to default (no restart)
2362!--    if the NAMELIST-parameter restart_time is omitted
2363       time_restart = 9999999.9
2364    ENDIF
2365
2366!
2367!-- Set default value of the time needed to terminate a model run
2368    IF ( termination_time_needed == -1.0 )  THEN
2369       IF ( host(1:3) == 'ibm' )  THEN
2370          termination_time_needed = 300.0
2371       ELSE
2372          termination_time_needed = 35.0
2373       ENDIF
2374    ENDIF
2375
2376!
2377!-- Check the time needed to terminate a model run
2378    IF ( host(1:3) == 't3e' )  THEN
2379!
2380!--    Time needed must be at least 30 seconds on all CRAY machines, because
2381!--    MPP_TREMAIN gives the remaining CPU time only in steps of 30 seconds
2382       IF ( termination_time_needed <= 30.0 )  THEN
2383          IF ( myid == 0 )  THEN
2384             PRINT*, '+++ check_parameters:  termination_time_needed', &
2385                      termination_time_needed
2386             PRINT*, '                       must be > 30.0 on host "', host, &
2387                     '"'
2388          ENDIF
2389          CALL local_stop
2390       ENDIF
2391    ELSEIF ( host(1:3) == 'ibm' )  THEN
2392!
2393!--    On IBM-regatta machines the time should be at least 300 seconds,
2394!--    because the job time consumed before executing palm (for compiling,
2395!--    copying of files, etc.) has to be regarded
2396       IF ( termination_time_needed < 300.0 )  THEN
2397          IF ( myid == 0 )  THEN
2398             PRINT*, '+++ WARNING: check_parameters:  termination_time_',  &
2399                         'needed', termination_time_needed
2400             PRINT*, '                                should be >= 300.0', &
2401                         ' on host "', host, '"'
2402          ENDIF
2403       ENDIF
2404    ENDIF
2405
2406
2407 END SUBROUTINE check_parameters
Note: See TracBrowser for help on using the repository browser.