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

Last change on this file since 1 was 1, checked in by raasch, 18 years ago

Initial repository layout and content

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