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

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

further changes concerning user-defined profiles

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