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

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

further preliminary uncomplete changes for ocean version

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