source: palm/tags/release-3.3/SOURCE/check_parameters.f90 @ 1983

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

updating comments and rc-file

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