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

Last change on this file since 141 was 90, checked in by raasch, 17 years ago

New:
---
Calculation and output of user-defined profiles. New &userpar parameters data_output_pr_user and max_pr_user.

check_parameters, flow_statistics, modules, parin, read_var_list, user_interface, write_var_list

Changed:


Division through dt_3d replaced by multiplication of the inverse. For performance optimisation, this is done in the loop calculating the divergence instead of using a seperate loop. (pres.f90) var_hom and var_sum renamed pr_palm.

data_output_profiles, flow_statistics, init_3d_model, modules, parin, pres, read_var_list, run_control, time_integration

Errors:


Bugfix: work_fft*_vec removed from some PRIVATE-declarations (poisfft).

Bugfix: field_chr renamed field_char (user_interface).

Bugfix: output of use_upstream_for_tke (header).

header, poisfft, user_interface

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