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

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

preliminary uncomplete changes for ocean version

  • Property svn:keywords set to Id
File size: 92.6 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 94 2007-06-01 15:25:22Z 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!-- In case of humidity or passive scalar, set boundary conditions for total
1008!-- water content / scalar
1009    IF ( humidity  .OR.  passive_scalar ) THEN
1010       IF ( humidity )  THEN
1011          sq = 'q'
1012       ELSE
1013          sq = 's'
1014       ENDIF
1015       IF ( bc_q_b == 'dirichlet' )  THEN
1016          ibc_q_b = 0
1017       ELSEIF ( bc_q_b == 'neumann' )  THEN
1018          ibc_q_b = 1
1019       ELSE
1020          IF ( myid == 0 )  THEN
1021             PRINT*, '+++ check_parameters:'
1022             PRINT*, '    unknown boundary condition: bc_', sq, '_b = ', bc_q_b
1023          ENDIF
1024          CALL local_stop
1025       ENDIF
1026       IF ( bc_q_t == 'dirichlet' )  THEN
1027          ibc_q_t = 0
1028       ELSEIF ( bc_q_t == 'neumann' )  THEN
1029          ibc_q_t = 1
1030       ELSE
1031          IF ( myid == 0 )  THEN
1032             PRINT*, '+++ check_parameters:'
1033             PRINT*, '    unknown boundary condition: bc_', sq, '_t = ', bc_q_t
1034          ENDIF
1035          CALL local_stop
1036       ENDIF
1037
1038       IF ( surface_waterflux == 0.0 )  constant_waterflux = .FALSE.
1039
1040!
1041!--    A given surface humidity implies Dirichlet boundary condition for
1042!--    humidity. In this case specification of a constant water flux is
1043!--    forbidden.
1044       IF ( ibc_q_b == 0  .AND.  constant_waterflux )  THEN
1045          IF ( myid == 0 )  THEN
1046             PRINT*, '+++ check_parameters:'
1047             PRINT*, '    boundary_condition: bc_', sq, '_b = ', bc_q_b
1048             PRINT*, '    is not allowed with prescribed surface flux'
1049          ENDIF
1050          CALL local_stop
1051       ENDIF
1052       IF ( constant_waterflux  .AND.  q_surface_initial_change /= 0.0 )  THEN
1053          IF ( myid == 0 )  THEN
1054             PRINT*, '+++ check_parameters: a prescribed surface flux is not'
1055             PRINT*, '    allowed with ', sq, '_surface_initial_change (/=0)', &
1056                     ' = ', q_surface_initial_change
1057          ENDIF
1058          CALL local_stop
1059       ENDIF
1060       
1061    ENDIF
1062
1063!
1064!-- Boundary conditions for horizontal components of wind speed
1065    IF ( bc_uv_b == 'dirichlet' )  THEN
1066       ibc_uv_b = 0
1067    ELSEIF ( bc_uv_b == 'neumann' )  THEN
1068       ibc_uv_b = 1
1069       IF ( prandtl_layer )  THEN
1070          IF ( myid == 0 )  THEN
1071             PRINT*, '+++ check_parameters:'
1072             PRINT*, '    boundary condition: bc_uv_b = ', TRIM( bc_uv_b ), &
1073                          ' is not allowed with'
1074             PRINT*, '    prandtl_layer = .TRUE.' 
1075          ENDIF
1076          CALL local_stop
1077       ENDIF
1078    ELSE
1079       IF ( myid == 0 )  THEN
1080          PRINT*, '+++ check_parameters:'
1081          PRINT*, '    unknown boundary condition: bc_uv_b = ', bc_uv_b
1082       ENDIF
1083       CALL local_stop
1084    ENDIF
1085    IF ( bc_uv_t == 'dirichlet' )  THEN
1086       ibc_uv_t = 0
1087    ELSEIF ( bc_uv_t == 'neumann' )  THEN
1088       ibc_uv_t = 1
1089    ELSE
1090       IF ( myid == 0 )  THEN
1091          PRINT*, '+++ check_parameters:'
1092          PRINT*, '    unknown boundary condition: bc_uv_t = ', bc_uv_t
1093       ENDIF
1094       CALL local_stop
1095    ENDIF
1096
1097!
1098!-- Compute and check, respectively, the Rayleigh Damping parameter
1099    IF ( rayleigh_damping_factor == -1.0 )  THEN
1100       IF ( momentum_advec == 'ups-scheme' )  THEN
1101          rayleigh_damping_factor = 0.01
1102       ELSE
1103          rayleigh_damping_factor = 0.0
1104       ENDIF
1105    ELSE
1106       IF ( rayleigh_damping_factor < 0.0 .OR. rayleigh_damping_factor > 1.0 ) &
1107       THEN
1108          IF ( myid == 0 )  THEN
1109             PRINT*, '+++ check_parameters:'
1110             PRINT*, '    rayleigh_damping_factor = ', rayleigh_damping_factor,&
1111                          ' out of range [0.0,1.0]'
1112          ENDIF
1113          CALL local_stop
1114       ENDIF
1115    ENDIF
1116
1117    IF ( rayleigh_damping_height == -1.0 )  THEN
1118       rayleigh_damping_height = 0.66666666666 * zu(nzt)
1119    ELSE
1120       IF ( rayleigh_damping_height < 0.0  .OR. &
1121            rayleigh_damping_height > zu(nzt) )  THEN
1122          IF ( myid == 0 )  THEN
1123             PRINT*, '+++ check_parameters:'
1124             PRINT*, '    rayleigh_damping_height = ', rayleigh_damping_height,&
1125                          ' out of range [0.0,', zu(nzt), ']'
1126          ENDIF
1127          CALL local_stop
1128       ENDIF
1129    ENDIF
1130
1131!
1132!-- Check limiters for Upstream-Spline scheme
1133    IF ( overshoot_limit_u < 0.0  .OR.  overshoot_limit_v < 0.0  .OR.  &
1134         overshoot_limit_w < 0.0  .OR.  overshoot_limit_pt < 0.0  .OR. &
1135         overshoot_limit_e < 0.0 )  THEN
1136       IF ( myid == 0 )  THEN
1137          PRINT*, '+++ check_parameters:'
1138          PRINT*, '    overshoot_limit_... < 0.0 is not allowed'
1139       ENDIF
1140       CALL local_stop
1141    ENDIF
1142    IF ( ups_limit_u < 0.0 .OR. ups_limit_v < 0.0 .OR. ups_limit_w < 0.0 .OR. &
1143         ups_limit_pt < 0.0 .OR. ups_limit_e < 0.0 )  THEN
1144       IF ( myid == 0 )  THEN
1145          PRINT*, '+++ check_parameters:'
1146          PRINT*, '    ups_limit_... < 0.0 is not allowed'
1147       ENDIF
1148       CALL local_stop
1149    ENDIF
1150
1151!
1152!-- Check number of chosen statistic regions. More than 10 regions are not
1153!-- allowed, because so far no more than 10 corresponding output files can
1154!-- be opened (cf. check_open)
1155    IF ( statistic_regions > 9  .OR.  statistic_regions < 0 )  THEN
1156       IF ( myid == 0 )  THEN
1157          PRINT*, '+++ check_parameters: Number of statistic_regions = ', &
1158                       statistic_regions+1
1159          PRINT*, '    Only 10 regions are allowed'
1160       ENDIF
1161       CALL local_stop
1162    ENDIF
1163    IF ( normalizing_region > statistic_regions  .OR. &
1164         normalizing_region < 0)  THEN
1165       IF ( myid == 0 )  THEN
1166          PRINT*, '+++ check_parameters: normalizing_region = ', &
1167                       normalizing_region, ' is unknown'
1168          PRINT*, '    Must be <= ', statistic_regions
1169       ENDIF
1170       CALL local_stop
1171    ENDIF
1172
1173!
1174!-- Set the default intervals for data output, if necessary
1175!-- NOTE: dt_dosp has already been set in package_parin
1176    IF ( dt_data_output /= 9999999.9 )  THEN
1177       IF ( dt_dopr           == 9999999.9 )  dt_dopr           = dt_data_output
1178       IF ( dt_dopts          == 9999999.9 )  dt_dopts          = dt_data_output
1179       IF ( dt_do2d_xy        == 9999999.9 )  dt_do2d_xy        = dt_data_output
1180       IF ( dt_do2d_xz        == 9999999.9 )  dt_do2d_xz        = dt_data_output
1181       IF ( dt_do2d_yz        == 9999999.9 )  dt_do2d_yz        = dt_data_output
1182       IF ( dt_do3d           == 9999999.9 )  dt_do3d           = dt_data_output
1183       IF ( dt_data_output_av == 9999999.9 )  dt_data_output_av = dt_data_output
1184    ENDIF
1185
1186!
1187!-- Set the default skip time intervals for data output, if necessary
1188    IF ( skip_time_dopr    == 9999999.9 ) &
1189                                       skip_time_dopr    = skip_time_data_output
1190    IF ( skip_time_dosp    == 9999999.9 ) &
1191                                       skip_time_dosp    = skip_time_data_output
1192    IF ( skip_time_do2d_xy == 9999999.9 ) &
1193                                       skip_time_do2d_xy = skip_time_data_output
1194    IF ( skip_time_do2d_xz == 9999999.9 ) &
1195                                       skip_time_do2d_xz = skip_time_data_output
1196    IF ( skip_time_do2d_yz == 9999999.9 ) &
1197                                       skip_time_do2d_yz = skip_time_data_output
1198    IF ( skip_time_do3d    == 9999999.9 ) &
1199                                       skip_time_do3d    = skip_time_data_output
1200    IF ( skip_time_data_output_av == 9999999.9 ) &
1201                                skip_time_data_output_av = skip_time_data_output
1202
1203!
1204!-- Check the average intervals (first for 3d-data, then for profiles and
1205!-- spectra)
1206    IF ( averaging_interval > dt_data_output_av )  THEN
1207       IF ( myid == 0 )  THEN
1208          PRINT*, '+++ check_parameters: average_interval=',              &
1209                       averaging_interval, ' must be <= dt_data_output=', &
1210                       dt_data_output
1211       ENDIF
1212       CALL local_stop
1213    ENDIF
1214
1215    IF ( averaging_interval_pr == 9999999.9 )  THEN
1216       averaging_interval_pr = averaging_interval
1217    ENDIF
1218
1219    IF ( averaging_interval_pr > dt_dopr )  THEN
1220       IF ( myid == 0 )  THEN
1221          PRINT*, '+++ check_parameters: averaging_interval_pr=', &
1222                       averaging_interval_pr, ' must be <= dt_dopr=', dt_dopr
1223       ENDIF
1224       CALL local_stop
1225    ENDIF
1226
1227    IF ( averaging_interval_sp == 9999999.9 )  THEN
1228       averaging_interval_sp = averaging_interval
1229    ENDIF
1230
1231    IF ( averaging_interval_sp > dt_dosp )  THEN
1232       IF ( myid == 0 )  THEN
1233          PRINT*, '+++ check_parameters: averaging_interval_sp=', &
1234                       averaging_interval_sp, ' must be <= dt_dosp=', dt_dosp
1235       ENDIF
1236       CALL local_stop
1237    ENDIF
1238
1239!
1240!-- Set the default interval for profiles entering the temporal average
1241    IF ( dt_averaging_input_pr == 9999999.9 )  THEN
1242       dt_averaging_input_pr = dt_averaging_input
1243    ENDIF
1244
1245!
1246!-- Set the default interval for the output of timeseries to a reasonable
1247!-- value (tries to minimize the number of calls of flow_statistics)
1248    IF ( dt_dots == 9999999.9 )  THEN
1249       IF ( averaging_interval_pr == 0.0 )  THEN
1250          dt_dots = MIN( dt_run_control, dt_dopr )
1251       ELSE
1252          dt_dots = MIN( dt_run_control, dt_averaging_input_pr )
1253       ENDIF
1254    ENDIF
1255
1256!
1257!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
1258    IF ( dt_averaging_input > averaging_interval )  THEN
1259       IF ( myid == 0 )  THEN
1260          PRINT*, '+++ check_parameters: dt_averaging_input=',                &
1261                       dt_averaging_input, ' must be <= averaging_interval=', &
1262                       averaging_interval
1263       ENDIF
1264       CALL local_stop
1265    ENDIF
1266
1267    IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
1268       IF ( myid == 0 )  THEN
1269          PRINT*, '+++ check_parameters: dt_averaging_input_pr=', &
1270                       dt_averaging_input_pr,                     &
1271                       ' must be <= averaging_interval_pr=',      &
1272                       averaging_interval_pr
1273       ENDIF
1274       CALL local_stop
1275    ENDIF
1276
1277!
1278!-- Set the default value for the integration interval of precipitation amount
1279    IF ( precipitation )  THEN
1280       IF ( precipitation_amount_interval == 9999999.9 )  THEN
1281          precipitation_amount_interval = dt_do2d_xy
1282       ELSE
1283          IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
1284             IF ( myid == 0 )  PRINT*, '+++ check_parameters: ',              &
1285                                       'precipitation_amount_interval =',     &
1286                                        precipitation_amount_interval,        &
1287                                       ' must not be larger than dt_do2d_xy', &
1288                                       ' = ', dt_do2d_xy   
1289       CALL local_stop
1290          ENDIF
1291       ENDIF
1292    ENDIF
1293
1294!
1295!-- Determine the number of output profiles and check whether they are
1296!-- permissible
1297    DO  WHILE ( data_output_pr(dopr_n+1) /= '          ' )
1298
1299       dopr_n = dopr_n + 1
1300       i = dopr_n
1301
1302!
1303!--    Determine internal profile number (for hom, homs)
1304!--    and store height levels
1305       SELECT CASE ( TRIM( data_output_pr(i) ) )
1306
1307          CASE ( 'u', '#u' )
1308             dopr_index(i) = 1
1309             dopr_unit(i)  = 'm/s'
1310             hom(:,2,1,:)  = SPREAD( zu, 2, statistic_regions+1 )
1311             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1312                dopr_initial_index(i) = 5
1313                hom(:,2,5,:)          = SPREAD( zu, 2, statistic_regions+1 )
1314                data_output_pr(i)     = data_output_pr(i)(2:)
1315             ENDIF
1316
1317          CASE ( 'v', '#v' )
1318             dopr_index(i) = 2
1319             dopr_unit(i)  = 'm/s'
1320             hom(:,2,2,:)  = SPREAD( zu, 2, statistic_regions+1 )
1321             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1322                dopr_initial_index(i) = 6
1323                hom(:,2,6,:)          = SPREAD( zu, 2, statistic_regions+1 )
1324                data_output_pr(i)     = data_output_pr(i)(2:)
1325             ENDIF
1326
1327          CASE ( 'w' )
1328             dopr_index(i) = 3
1329             dopr_unit(i)  = 'm/s'
1330             hom(:,2,3,:)  = SPREAD( zw, 2, statistic_regions+1 )
1331
1332          CASE ( 'pt', '#pt' )
1333             IF ( .NOT. cloud_physics ) THEN
1334                dopr_index(i) = 4
1335                dopr_unit(i)  = 'K'
1336                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
1337                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1338                   dopr_initial_index(i) = 7
1339                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
1340                   hom(nzb,2,7,:)        = 0.0    ! because zu(nzb) is negative
1341                   data_output_pr(i)     = data_output_pr(i)(2:)
1342                ENDIF
1343             ELSE
1344                dopr_index(i) = 43
1345                dopr_unit(i)  = 'K'
1346                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
1347                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1348                   dopr_initial_index(i) = 28
1349                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
1350                   hom(nzb,2,28,:)       = 0.0    ! because zu(nzb) is negative
1351                   data_output_pr(i)     = data_output_pr(i)(2:)
1352                ENDIF
1353             ENDIF
1354
1355          CASE ( 'e' )
1356             dopr_index(i)  = 8
1357             dopr_unit(i)   = 'm2/s2'
1358             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
1359             hom(nzb,2,8,:) = 0.0
1360
1361          CASE ( 'km', '#km' )
1362             dopr_index(i)  = 9
1363             dopr_unit(i)   = 'm2/s'
1364             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
1365             hom(nzb,2,9,:) = 0.0
1366             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1367                dopr_initial_index(i) = 23
1368                hom(:,2,23,:)         = hom(:,2,9,:)
1369                data_output_pr(i)     = data_output_pr(i)(2:)
1370             ENDIF
1371
1372          CASE ( 'kh', '#kh' )
1373             dopr_index(i)   = 10
1374             dopr_unit(i)    = 'm2/s'
1375             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
1376             hom(nzb,2,10,:) = 0.0
1377             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1378                dopr_initial_index(i) = 24
1379                hom(:,2,24,:)         = hom(:,2,10,:)
1380                data_output_pr(i)     = data_output_pr(i)(2:)
1381             ENDIF
1382
1383          CASE ( 'l', '#l' )
1384             dopr_index(i)   = 11
1385             dopr_unit(i)    = 'm'
1386             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
1387             hom(nzb,2,11,:) = 0.0
1388             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1389                dopr_initial_index(i) = 25
1390                hom(:,2,25,:)         = hom(:,2,11,:)
1391                data_output_pr(i)     = data_output_pr(i)(2:)
1392             ENDIF
1393
1394          CASE ( 'w"u"' )
1395             dopr_index(i) = 12
1396             dopr_unit(i)  = 'm2/s2'
1397             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
1398             IF ( prandtl_layer )  hom(nzb,2,12,:) = zu(1)
1399
1400          CASE ( 'w*u*' )
1401             dopr_index(i) = 13
1402             dopr_unit(i)  = 'm2/s2'
1403             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
1404
1405          CASE ( 'w"v"' )
1406             dopr_index(i) = 14
1407             dopr_unit(i)  = 'm2/s2'
1408             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
1409             IF ( prandtl_layer )  hom(nzb,2,14,:) = zu(1)
1410
1411          CASE ( 'w*v*' )
1412             dopr_index(i) = 15
1413             dopr_unit(i)  = 'm2/s2'
1414             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
1415
1416          CASE ( 'w"pt"' )
1417             dopr_index(i) = 16
1418             dopr_unit(i)  = 'K m/s'
1419             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
1420
1421          CASE ( 'w*pt*' )
1422             dopr_index(i) = 17
1423             dopr_unit(i)  = 'K m/s'
1424             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
1425
1426          CASE ( 'wpt' )
1427             dopr_index(i) = 18
1428             dopr_unit(i)  = 'K m/s'
1429             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
1430
1431          CASE ( 'wu' )
1432             dopr_index(i) = 19
1433             dopr_unit(i)  = 'm2/s2'
1434             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
1435             IF ( prandtl_layer )  hom(nzb,2,19,:) = zu(1)
1436
1437          CASE ( 'wv' )
1438             dopr_index(i) = 20
1439             dopr_unit(i)  = 'm2/s2'
1440             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
1441             IF ( prandtl_layer )  hom(nzb,2,20,:) = zu(1)
1442
1443          CASE ( 'w*pt*BC' )
1444             dopr_index(i) = 21
1445             dopr_unit(i)  = 'K m/s'
1446             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
1447
1448          CASE ( 'wptBC' )
1449             dopr_index(i) = 22
1450             dopr_unit(i)  = 'K m/s'
1451             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
1452
1453          CASE ( 'u*2' )
1454             dopr_index(i) = 30
1455             dopr_unit(i)  = 'm2/s2'
1456             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
1457
1458          CASE ( 'v*2' )
1459             dopr_index(i) = 31
1460             dopr_unit(i)  = 'm2/s2'
1461             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
1462
1463          CASE ( 'w*2' )
1464             dopr_index(i) = 32
1465             dopr_unit(i)  = 'm2/s2'
1466             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
1467
1468          CASE ( 'pt*2' )
1469             dopr_index(i) = 33
1470             dopr_unit(i)  = 'K2'
1471             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
1472
1473          CASE ( 'e*' )
1474             dopr_index(i) = 34
1475             dopr_unit(i)  = 'm2/s2'
1476             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
1477
1478          CASE ( 'w*2pt*' )
1479             dopr_index(i) = 35
1480             dopr_unit(i)  = 'K m2/s2'
1481             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
1482
1483          CASE ( 'w*pt*2' )
1484             dopr_index(i) = 36
1485             dopr_unit(i)  = 'K2 m/s'
1486             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
1487
1488          CASE ( 'w*e*' )
1489             dopr_index(i) = 37
1490             dopr_unit(i)  = 'm3/s3'
1491             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
1492
1493          CASE ( 'w*3' )
1494             dopr_index(i) = 38
1495             dopr_unit(i)  = 'm3/s3'
1496             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
1497
1498          CASE ( 'Sw' )
1499             dopr_index(i) = 39
1500             dopr_unit(i)  = 'none'
1501             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
1502
1503          CASE ( 'q', '#q' )
1504             IF ( .NOT. cloud_physics )  THEN
1505                IF ( myid == 0 )  THEN
1506                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1507                           data_output_pr(i),                          &
1508                           '    is not implemented for cloud_physics = FALSE'
1509                ENDIF
1510                CALL local_stop
1511             ELSE
1512                dopr_index(i) = 41
1513                dopr_unit(i)  = 'kg/kg'
1514                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
1515                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1516                   dopr_initial_index(i) = 26
1517                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
1518                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
1519                   data_output_pr(i)     = data_output_pr(i)(2:)
1520                ENDIF
1521             ENDIF
1522
1523          CASE ( 's', '#s' )
1524             IF ( .NOT. passive_scalar )  THEN
1525                IF ( myid == 0 )  THEN
1526                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1527                           data_output_pr(i),                          &
1528                           '    is not implemented for passive_scalar = FALSE'
1529                ENDIF
1530                CALL local_stop
1531             ELSE
1532                dopr_index(i) = 41
1533                dopr_unit(i)  = 'kg/m3'
1534                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
1535                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1536                   dopr_initial_index(i) = 26
1537                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
1538                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
1539                   data_output_pr(i)     = data_output_pr(i)(2:)
1540                ENDIF
1541             ENDIF
1542
1543          CASE ( 'qv', '#qv' )
1544             IF ( .NOT. cloud_physics ) THEN
1545                dopr_index(i) = 41
1546                dopr_unit(i)  = 'kg/kg'
1547                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
1548                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1549                   dopr_initial_index(i) = 26
1550                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
1551                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
1552                   data_output_pr(i)     = data_output_pr(i)(2:)
1553                ENDIF
1554             ELSE
1555                dopr_index(i) = 42
1556                dopr_unit(i)  = 'kg/kg'
1557                hom(:,2,42,:) = SPREAD( zu, 2, statistic_regions+1 )
1558                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1559                   dopr_initial_index(i) = 27
1560                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
1561                   hom(nzb,2,27,:)       = 0.0    ! weil zu(nzb) negativ ist
1562                   data_output_pr(i)     = data_output_pr(i)(2:)
1563                ENDIF
1564             ENDIF
1565
1566          CASE ( 'lpt', '#lpt' )
1567             IF ( .NOT. cloud_physics ) THEN
1568                IF ( myid == 0 )  THEN
1569                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1570                           data_output_pr(i),                          &
1571                           '    is not implemented for cloud_physics = FALSE'
1572                ENDIF
1573                CALL local_stop
1574             ELSE
1575                dopr_index(i) = 4
1576                dopr_unit(i)  = 'K'
1577                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
1578                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1579                   dopr_initial_index(i) = 7
1580                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
1581                   hom(nzb,2,7,:)        = 0.0    ! weil zu(nzb) negativ ist
1582                   data_output_pr(i)     = data_output_pr(i)(2:)
1583                ENDIF
1584             ENDIF
1585
1586          CASE ( 'vpt', '#vpt' )
1587             dopr_index(i) = 44
1588             dopr_unit(i)  = 'K'
1589             hom(:,2,44,:) = SPREAD( zu, 2, statistic_regions+1 )
1590             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1591                dopr_initial_index(i) = 29
1592                hom(:,2,29,:)         = SPREAD( zu, 2, statistic_regions+1 )
1593                hom(nzb,2,29,:)       = 0.0    ! weil zu(nzb) negativ ist
1594                data_output_pr(i)     = data_output_pr(i)(2:)
1595             ENDIF
1596
1597          CASE ( 'w"vpt"' )
1598             dopr_index(i) = 45
1599             dopr_unit(i)  = 'K m/s'
1600             hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
1601
1602          CASE ( 'w*vpt*' )
1603             dopr_index(i) = 46
1604             dopr_unit(i)  = 'K m/s'
1605             hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
1606
1607          CASE ( 'wvpt' )
1608             dopr_index(i) = 47
1609             dopr_unit(i)  = 'K m/s'
1610             hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
1611
1612          CASE ( 'w"q"' )
1613             IF ( .NOT. cloud_physics ) THEN
1614                IF ( myid == 0 )  THEN
1615                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1616                           data_output_pr(i),                          &
1617                           '    is not implemented for cloud_physics = FALSE'
1618                ENDIF
1619                CALL local_stop
1620             ELSE
1621                dopr_index(i) = 48
1622                dopr_unit(i)  = 'kg/kg m/s'
1623                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
1624             ENDIF
1625
1626          CASE ( 'w*q*' )
1627             IF ( .NOT. cloud_physics ) THEN
1628                IF ( myid == 0 )  THEN
1629                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1630                           data_output_pr(i),                          &
1631                           '    is not implemented for cloud_physics = FALSE'
1632                ENDIF
1633                CALL local_stop
1634             ELSE
1635                dopr_index(i) = 49
1636                dopr_unit(i)  = 'kg/kg m/s'
1637                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
1638             ENDIF
1639
1640          CASE ( 'wq' )
1641             IF ( .NOT. cloud_physics ) THEN
1642                IF ( myid == 0 )  THEN
1643                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1644                           data_output_pr(i),                          &
1645                           '    is not implemented for cloud_physics = FALSE'
1646                ENDIF
1647                CALL local_stop
1648             ELSE
1649                dopr_index(i) = 50
1650                dopr_unit(i)  = 'kg/kg m/s'
1651                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
1652             ENDIF
1653
1654          CASE ( 'w"s"' )
1655             IF ( .NOT. passive_scalar ) THEN
1656                IF ( myid == 0 )  THEN
1657                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1658                           data_output_pr(i),                          &
1659                           '    is not implemented for passive_scalar = FALSE'
1660                ENDIF
1661                CALL local_stop
1662             ELSE
1663                dopr_index(i) = 48
1664                dopr_unit(i)  = 'kg/m3 m/s'
1665                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
1666             ENDIF
1667
1668          CASE ( 'w*s*' )
1669             IF ( .NOT. passive_scalar ) THEN
1670                IF ( myid == 0 )  THEN
1671                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1672                           data_output_pr(i),                          &
1673                           '    is not implemented for passive_scalar = FALSE'
1674                ENDIF
1675                CALL local_stop
1676             ELSE
1677                dopr_index(i) = 49
1678                dopr_unit(i)  = 'kg/m3 m/s'
1679                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
1680             ENDIF
1681
1682          CASE ( 'ws' )
1683             IF ( .NOT. passive_scalar ) THEN
1684                IF ( myid == 0 )  THEN
1685                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1686                           data_output_pr(i),                          &
1687                           '    is not implemented for passive_scalar = FALSE'
1688                ENDIF
1689                CALL local_stop
1690             ELSE
1691                dopr_index(i) = 50
1692                dopr_unit(i)  = 'kg/m3 m/s'
1693                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
1694             ENDIF
1695
1696          CASE ( 'w"qv"' )
1697             IF ( humidity  .AND.  .NOT. cloud_physics ) &
1698             THEN
1699                dopr_index(i) = 48
1700                dopr_unit(i)  = 'kg/kg m/s'
1701                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
1702             ELSEIF( humidity .AND. cloud_physics ) THEN
1703                dopr_index(i) = 51
1704                dopr_unit(i)  = 'kg/kg m/s'
1705                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
1706             ELSE
1707                IF ( myid == 0 )  THEN
1708                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1709                           data_output_pr(i),                          &
1710                           '    is not implemented for cloud_physics = FALSE', &
1711                           '    and                    humidity      = FALSE'
1712                ENDIF
1713                CALL local_stop                   
1714             ENDIF
1715
1716          CASE ( 'w*qv*' )
1717             IF ( humidity  .AND.  .NOT. cloud_physics ) &
1718             THEN
1719                dopr_index(i) = 49
1720                dopr_unit(i)  = 'kg/kg m/s'
1721                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
1722             ELSEIF( humidity .AND. cloud_physics ) THEN
1723                dopr_index(i) = 52
1724                dopr_unit(i)  = 'kg/kg m/s'
1725                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
1726             ELSE
1727                IF ( myid == 0 )  THEN
1728                   PRINT*, '+++ check_parameters:  data_output_pr = ',         &
1729                           data_output_pr(i),                                  &
1730                           '    is not implemented for cloud_physics = FALSE', &
1731                           '                       and humidity      = FALSE'
1732                ENDIF
1733                CALL local_stop                   
1734             ENDIF
1735
1736          CASE ( 'wqv' )
1737             IF ( humidity  .AND.  .NOT. cloud_physics ) &
1738             THEN
1739                dopr_index(i) = 50
1740                dopr_unit(i)  = 'kg/kg m/s'
1741                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
1742             ELSEIF( humidity .AND. cloud_physics ) THEN
1743                dopr_index(i) = 53
1744                dopr_unit(i)  = 'kg/kg m/s'
1745                hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
1746             ELSE
1747                IF ( myid == 0 )  THEN
1748                   PRINT*, '+++ check_parameters:  data_output_pr = ',         &
1749                           data_output_pr(i),                                  &
1750                           '    is not implemented for cloud_physics = FALSE', &
1751                           '                       and humidity      = FALSE'
1752                ENDIF
1753                CALL local_stop                   
1754             ENDIF
1755
1756          CASE ( 'ql' )
1757             IF ( .NOT. cloud_physics  .AND.  .NOT. cloud_droplets )  THEN
1758                IF ( myid == 0 )  THEN
1759                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
1760                           data_output_pr(i),                          &
1761                           '    is not implemented for cloud_physics = FALSE'
1762                ENDIF
1763                CALL local_stop
1764             ELSE
1765                dopr_index(i) = 54
1766                dopr_unit(i)  = 'kg/kg'
1767                hom(:,2,54,:)  = SPREAD( zu, 2, statistic_regions+1 )
1768             ENDIF
1769
1770          CASE ( 'w*u*u*/dz' )
1771             dopr_index(i) = 55
1772             dopr_unit(i)  = 'm2/s3'
1773             hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
1774
1775          CASE ( 'w*p*/dz' )
1776             dopr_index(i) = 56
1777             dopr_unit(i)  = 'm2/s3'
1778             hom(:,2,56,:) = SPREAD( zu, 2, statistic_regions+1 )
1779
1780          CASE ( 'w"e/dz' )
1781             dopr_index(i) = 57
1782             dopr_unit(i)  = 'm2/s3'
1783             hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
1784
1785          CASE ( 'u"pt"' )
1786             dopr_index(i) = 58
1787             dopr_unit(i)  = 'K m/s'
1788             hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
1789
1790          CASE ( 'u*pt*' )
1791             dopr_index(i) = 59
1792             dopr_unit(i)  = 'K m/s'
1793             hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
1794
1795          CASE ( 'upt_t' )
1796             dopr_index(i) = 60
1797             dopr_unit(i)  = 'K m/s'
1798             hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
1799
1800          CASE ( 'v"pt"' )
1801             dopr_index(i) = 61
1802             dopr_unit(i)  = 'K m/s'
1803             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
1804             
1805          CASE ( 'v*pt*' )
1806             dopr_index(i) = 62
1807             dopr_unit(i)  = 'K m/s'
1808             hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
1809
1810          CASE ( 'vpt_t' )
1811             dopr_index(i) = 63
1812             dopr_unit(i)  = 'K m/s'
1813             hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
1814
1815
1816          CASE DEFAULT
1817
1818             CALL user_check_data_output_pr( data_output_pr(i), i, unit )
1819
1820             IF ( unit == 'illegal' )  THEN
1821                IF ( myid == 0 )  THEN
1822                   IF ( data_output_pr_user(1) /= ' ' )  THEN
1823                      PRINT*, '+++ check_parameters:  illegal value for data_',&
1824                                   'output_pr or data_output_pr_user: "',      &
1825                                   TRIM( data_output_pr(i) ), '"'
1826                   ELSE
1827                      PRINT*, '+++ check_parameters:  illegal value for data_',&
1828                                   'output_pr: "', TRIM( data_output_pr(i) ),'"'
1829                   ENDIF
1830                ENDIF
1831                CALL local_stop
1832             ENDIF
1833
1834       END SELECT
1835!
1836!--    Check to which of the predefined coordinate systems the profile belongs
1837       DO  k = 1, crmax
1838          IF ( INDEX( cross_profiles(k), ' '//TRIM( data_output_pr(i) )//' ' ) &
1839               /=0 ) &
1840          THEN
1841             dopr_crossindex(i) = k
1842             EXIT
1843          ENDIF
1844       ENDDO
1845!
1846!--    Generate the text for the labels of the PROFIL output file. "-characters
1847!--    must be substituted, otherwise PROFIL would interpret them as TeX
1848!--    control characters
1849       dopr_label(i) = data_output_pr(i)
1850       position = INDEX( dopr_label(i) , '"' )
1851       DO WHILE ( position /= 0 )
1852          dopr_label(i)(position:position) = ''''
1853          position = INDEX( dopr_label(i) , '"' )
1854       ENDDO
1855
1856    ENDDO
1857
1858!
1859!-- y-value range of the coordinate system (PROFIL).
1860!-- x-value range determined in plot_1d.
1861    IF ( .NOT. ocean )  THEN
1862       cross_uymin = 0.0
1863       IF ( z_max_do1d == -1.0 )  THEN
1864          cross_uymax = zu(nzt+1)
1865       ELSEIF ( z_max_do1d < zu(nzb+1)  .OR.  z_max_do1d > zu(nzt+1) )  THEN
1866          IF ( myid == 0 )  PRINT*, '+++ check_parameters:  z_max_do1d=',  &
1867                                    z_max_do1d,' must be >= ', zu(nzb+1),  &
1868                                    ' or <= ', zu(nzt+1)
1869          CALL local_stop
1870       ELSE
1871          cross_uymax = z_max_do1d
1872       ENDIF
1873    ENDIF
1874
1875!
1876!-- Check whether the chosen normalizing factor for the coordinate systems is
1877!-- permissible
1878    DO  i = 1, crmax
1879       SELECT CASE ( TRIM( cross_normalized_x(i) ) )  ! TRIM required on IBM
1880
1881          CASE ( '', 'wpt0', 'ws2', 'tsw2', 'ws3', 'ws2tsw', 'wstsw2' )
1882             j = 0
1883
1884          CASE DEFAULT
1885             IF ( myid == 0 )  THEN
1886                PRINT*, '+++ check_parameters: unknown normalize method'
1887                PRINT*, '    cross_normalized_x="',cross_normalized_x(i),'"'
1888             ENDIF
1889             CALL local_stop
1890
1891       END SELECT
1892       SELECT CASE ( TRIM( cross_normalized_y(i) ) )  ! TRIM required on IBM
1893
1894          CASE ( '', 'z_i' )
1895             j = 0
1896
1897          CASE DEFAULT
1898             IF ( myid == 0 )  THEN
1899                PRINT*, '+++ check_parameters: unknown normalize method'
1900                PRINT*, '    cross_normalized_y="',cross_normalized_y(i),'"'
1901             ENDIF
1902             CALL local_stop
1903
1904       END SELECT
1905    ENDDO
1906!
1907!-- Check normalized y-value range of the coordinate system (PROFIL)
1908    IF ( z_max_do1d_normalized /= -1.0  .AND.  z_max_do1d_normalized <= 0.0 ) &
1909    THEN
1910       IF ( myid == 0 )  PRINT*,'+++ check_parameters:  z_max_do1d_normalize', &
1911                                'd=', z_max_do1d_normalized, ' must be >= 0.0'
1912       CALL local_stop
1913    ENDIF
1914
1915
1916!
1917!-- Append user-defined data output variables to the standard data output
1918    IF ( data_output_user(1) /= ' ' )  THEN
1919       i = 1
1920       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
1921          i = i + 1
1922       ENDDO
1923       j = 1
1924       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 100 )
1925          IF ( i > 100 )  THEN
1926             IF ( myid == 0 )  THEN
1927                PRINT*, '+++ check_parameters: number of output quantitities', &
1928                             ' given by data_output and data_output_user'
1929                PRINT*, '                      exceeds the limit of 100'
1930             ENDIF
1931             CALL local_stop
1932          ENDIF
1933          data_output(i) = data_output_user(j)
1934          i = i + 1
1935          j = j + 1
1936       ENDDO
1937    ENDIF
1938
1939!
1940!-- Check and set steering parameters for 2d/3d data output and averaging
1941    i   = 1
1942    DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
1943!
1944!--    Check for data averaging
1945       ilen = LEN_TRIM( data_output(i) )
1946       j = 0                                                 ! no data averaging
1947       IF ( ilen > 3 )  THEN
1948          IF ( data_output(i)(ilen-2:ilen) == '_av' )  THEN
1949             j = 1                                           ! data averaging
1950             data_output(i) = data_output(i)(1:ilen-3)
1951          ENDIF
1952       ENDIF
1953!
1954!--    Check for cross section or volume data
1955       ilen = LEN_TRIM( data_output(i) )
1956       k = 0                                                   ! 3d data
1957       var = data_output(i)(1:ilen)
1958       IF ( ilen > 3 )  THEN
1959          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR. &
1960               data_output(i)(ilen-2:ilen) == '_xz'  .OR. &
1961               data_output(i)(ilen-2:ilen) == '_yz' )  THEN
1962             k = 1                                             ! 2d data
1963             var = data_output(i)(1:ilen-3)
1964          ENDIF
1965       ENDIF
1966!
1967!--    Check for allowed value and set units
1968       SELECT CASE ( TRIM( var ) )
1969
1970          CASE ( 'e' )
1971             IF ( constant_diffusion )  THEN
1972                IF ( myid == 0 )  THEN
1973                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1974                                '" requires constant_diffusion = .FALSE.'
1975                ENDIF
1976                CALL local_stop
1977             ENDIF
1978             unit = 'm2/s2'
1979
1980          CASE ( 'pc', 'pr' )
1981             IF ( .NOT. particle_advection )  THEN
1982                IF ( myid == 0 )  THEN
1983                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1984                                '" requires particle package'
1985                   PRINT*, '                      (mrun-option "-p particles")'
1986                ENDIF
1987                CALL local_stop
1988             ENDIF
1989             IF ( TRIM( var ) == 'pc' )  unit = 'number'
1990             IF ( TRIM( var ) == 'pr' )  unit = 'm'
1991
1992          CASE ( 'q', 'vpt' )
1993             IF ( .NOT. humidity )  THEN
1994                IF ( myid == 0 )  THEN
1995                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
1996                                '" requires humidity = .TRUE.'
1997                ENDIF
1998                CALL local_stop
1999             ENDIF
2000             IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
2001             IF ( TRIM( var ) == 'vpt' )  unit = 'K'
2002
2003          CASE ( 'ql' )
2004             IF ( .NOT. ( cloud_physics  .OR.  cloud_droplets ) )  THEN
2005                IF ( myid == 0 )  THEN
2006                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
2007                                '" requires cloud_physics = .TRUE.'
2008                   PRINT*, '                      or cloud_droplets = .TRUE.'
2009                ENDIF
2010                CALL local_stop
2011             ENDIF
2012             unit = 'kg/kg'
2013
2014          CASE ( 'ql_c', 'ql_v', 'ql_vp' )
2015             IF ( .NOT. cloud_droplets )  THEN
2016                IF ( myid == 0 )  THEN
2017                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
2018                                '" requires cloud_droplets = .TRUE.'
2019                ENDIF
2020                CALL local_stop
2021             ENDIF
2022             IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
2023             IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
2024             IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
2025
2026          CASE ( 'qv' )
2027             IF ( .NOT. cloud_physics )  THEN
2028                IF ( myid == 0 )  THEN
2029                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
2030                                '" requires cloud_physics = .TRUE.'
2031                ENDIF
2032                CALL local_stop
2033             ENDIF
2034             unit = 'kg/kg'
2035
2036          CASE ( 's' )
2037             IF ( .NOT. passive_scalar )  THEN
2038                IF ( myid == 0 )  THEN
2039                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
2040                                '" requires passive_scalar = .TRUE.'
2041                ENDIF
2042                CALL local_stop
2043             ENDIF
2044             unit = 'conc'
2045
2046          CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*' )
2047             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
2048                IF ( myid == 0 )  THEN
2049                   PRINT*, '+++ check_parameters:  illegal value for data_',&
2050                                'output: "', TRIM( var ), '" is only allowed'
2051                   PRINT*, '                       for horizontal cross section'
2052                ENDIF
2053                CALL local_stop
2054             ENDIF
2055             IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
2056                IF ( myid == 0 )  THEN
2057                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
2058                                '" requires cloud_physics = .TRUE.'
2059                ENDIF
2060                CALL local_stop
2061             ENDIF
2062             IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
2063                IF ( myid == 0 )  THEN
2064                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
2065                                '" requires precipitation = .TRUE.'
2066                ENDIF
2067                CALL local_stop
2068             ENDIF
2069             IF ( TRIM( var ) == 'pra*'  .AND.  j == 1 )  THEN
2070                IF ( myid == 0 )  THEN
2071                   PRINT*, '+++ check_parameters: temporal averaging of ', &
2072                           ' precipitation amount "', TRIM( var ),         &
2073                           '" not possible' 
2074                ENDIF
2075                CALL local_stop
2076             ENDIF
2077             IF ( TRIM( var ) == 'prr*'  .AND.  .NOT. precipitation )  THEN
2078                IF ( myid == 0 )  THEN
2079                   PRINT*, '+++ check_parameters: output of "', TRIM( var ), &
2080                                '" requires precipitation = .TRUE.'
2081                ENDIF
2082                CALL local_stop
2083             ENDIF
2084
2085
2086             IF ( TRIM( var ) == 'u*'   )  unit = 'm/s'
2087             IF ( TRIM( var ) == 't*'   )  unit = 'K'
2088             IF ( TRIM( var ) == 'lwp*' )  unit = 'kg/kg*m'
2089             IF ( TRIM( var ) == 'pra*' )  unit = 'mm'
2090             IF ( TRIM( var ) == 'prr*' )  unit = 'mm/s'
2091             IF ( TRIM( var ) == 'z0*'  )  unit = 'm'
2092
2093          CASE ( 'p', 'pt', 'u', 'v', 'w' )
2094             IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
2095             IF ( TRIM( var ) == 'pt' )  unit = 'K'
2096             IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
2097             IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
2098             IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
2099             CONTINUE
2100
2101          CASE DEFAULT
2102             CALL user_check_data_output( var, unit )
2103
2104             IF ( unit == 'illegal' )  THEN
2105                IF ( myid == 0 )  THEN
2106                   IF ( data_output_user(1) /= ' ' )  THEN
2107                      PRINT*, '+++ check_parameters:  illegal value for data_',&
2108                                   'output or data_output_user: "',            &
2109                                   TRIM( data_output(i) ), '"'
2110                   ELSE
2111                      PRINT*, '+++ check_parameters:  illegal value for data_',&
2112                                   'output: "', TRIM( data_output(i) ), '"'
2113                   ENDIF
2114                ENDIF
2115                CALL local_stop
2116             ENDIF
2117
2118       END SELECT
2119!
2120!--    Set the internal steering parameters appropriately
2121       IF ( k == 0 )  THEN
2122          do3d_no(j)              = do3d_no(j) + 1
2123          do3d(j,do3d_no(j))      = data_output(i)
2124          do3d_unit(j,do3d_no(j)) = unit
2125       ELSE
2126          do2d_no(j)              = do2d_no(j) + 1
2127          do2d(j,do2d_no(j))      = data_output(i)
2128          do2d_unit(j,do2d_no(j)) = unit
2129          IF ( data_output(i)(ilen-2:ilen) == '_xy' )  THEN
2130             data_output_xy(j) = .TRUE.
2131          ENDIF
2132          IF ( data_output(i)(ilen-2:ilen) == '_xz' )  THEN
2133             data_output_xz(j) = .TRUE.
2134          ENDIF
2135          IF ( data_output(i)(ilen-2:ilen) == '_yz' )  THEN
2136             data_output_yz(j) = .TRUE.
2137          ENDIF
2138       ENDIF
2139
2140       IF ( j == 1 )  THEN
2141!
2142!--       Check, if variable is already subject to averaging
2143          found = .FALSE.
2144          DO  k = 1, doav_n
2145             IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
2146          ENDDO
2147
2148          IF ( .NOT. found )  THEN
2149             doav_n = doav_n + 1
2150             doav(doav_n) = var
2151          ENDIF
2152       ENDIF
2153
2154       i = i + 1
2155    ENDDO
2156
2157!
2158!-- Store sectional planes in one shared array
2159    section(:,1) = section_xy
2160    section(:,2) = section_xz
2161    section(:,3) = section_yz
2162
2163!
2164!-- Upper plot limit (grid point value) for 1D profiles
2165    IF ( z_max_do1d == -1.0 )  THEN
2166       nz_do1d = nzt+1
2167    ELSE
2168       DO  k = nzb+1, nzt+1
2169          nz_do1d = k
2170          IF ( zw(k) > z_max_do1d )  EXIT
2171       ENDDO
2172    ENDIF
2173
2174!
2175!-- Upper plot limit for 2D vertical sections
2176    IF ( z_max_do2d == -1.0 )  z_max_do2d = zu(nzt)
2177    IF ( z_max_do2d < zu(nzb+1)  .OR.  z_max_do2d > zu(nzt) )  THEN
2178       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  z_max_do2d=',        &
2179                                 z_max_do2d, ' must be >= ', zu(nzb+1),       &
2180                                 '(zu(nzb+1)) and <= ', zu(nzt), ' (zu(nzt))'
2181       CALL local_stop
2182    ENDIF
2183
2184!
2185!-- Upper plot limit for 3D arrays
2186    IF ( nz_do3d == -9999 )  nz_do3d = nzt + 1
2187
2188!
2189!-- Determine and check accuracy for compressed 3D plot output
2190    IF ( do3d_compress )  THEN
2191!
2192!--    Compression only permissible on T3E machines
2193       IF ( host(1:3) /= 't3e' )  THEN
2194          IF ( myid == 0 )  THEN
2195             PRINT*, '+++ check_parameters: do3d_compress = .TRUE. not allow', &
2196                          'ed on host "', TRIM( host ), '"'
2197          ENDIF
2198          CALL local_stop
2199       ENDIF
2200
2201       i = 1
2202       DO  WHILE ( do3d_comp_prec(i) /= ' ' )
2203
2204          ilen = LEN_TRIM( do3d_comp_prec(i) )
2205          IF ( LLT( do3d_comp_prec(i)(ilen:ilen), '0' ) .OR. &
2206               LGT( do3d_comp_prec(i)(ilen:ilen), '9' ) )  THEN
2207             IF ( myid == 0 )  THEN
2208                PRINT*, '+++ check_parameters: illegal precision: ', &
2209                        'do3d_comp_prec(', i, ')="', TRIM(do3d_comp_prec(i)),'"'
2210             ENDIF
2211             CALL local_stop
2212          ENDIF
2213
2214          prec = IACHAR( do3d_comp_prec(i)(ilen:ilen) ) - IACHAR( '0' )
2215          var = do3d_comp_prec(i)(1:ilen-1)
2216
2217          SELECT CASE ( var )
2218
2219             CASE ( 'u' )
2220                j = 1
2221             CASE ( 'v' )
2222                j = 2
2223             CASE ( 'w' )
2224                j = 3
2225             CASE ( 'p' )
2226                j = 4
2227             CASE ( 'pt' )
2228                j = 5
2229
2230             CASE DEFAULT
2231                IF ( myid == 0 )  THEN
2232                   PRINT*, '+++ check_parameters: unknown variable in ', &
2233                               'assignment'
2234                   PRINT*, '    do3d_comp_prec(', i, ')="', &
2235                           TRIM( do3d_comp_prec(i) ),'"'
2236                ENDIF
2237                CALL local_stop               
2238
2239          END SELECT
2240
2241          plot_3d_precision(j)%precision = prec
2242          i = i + 1
2243
2244       ENDDO
2245    ENDIF
2246
2247!
2248!-- Check the data output format(s)
2249    IF ( data_output_format(1) == ' ' )  THEN
2250!
2251!--    Default value
2252       netcdf_output = .TRUE.
2253    ELSE
2254       i = 1
2255       DO  WHILE ( data_output_format(i) /= ' ' )
2256
2257          SELECT CASE ( data_output_format(i) )
2258
2259             CASE ( 'netcdf' )
2260                netcdf_output = .TRUE.
2261             CASE ( 'iso2d' )
2262                iso2d_output  = .TRUE.
2263             CASE ( 'profil' )
2264                profil_output = .TRUE.
2265             CASE ( 'avs' )
2266                avs_output    = .TRUE.
2267
2268             CASE DEFAULT
2269                IF ( myid == 0 )  THEN
2270                   PRINT*, '+++ check_parameters:'
2271                   PRINT*, '    unknown value for data_output_format "', &
2272                                TRIM( data_output_format(i) ),'"'
2273                ENDIF
2274                CALL local_stop               
2275
2276          END SELECT
2277
2278          i = i + 1
2279          IF ( i > 10 )  EXIT
2280
2281       ENDDO
2282
2283    ENDIF
2284
2285!
2286!-- Check netcdf precison
2287    ldum = .FALSE.
2288    CALL define_netcdf_header( 'ch', ldum, 0 )
2289
2290!
2291!-- Check, whether a constant diffusion coefficient shall be used
2292    IF ( km_constant /= -1.0 )  THEN
2293       IF ( km_constant < 0.0 )  THEN
2294          IF ( myid == 0 )  PRINT*, '+++ check_parameters:  km_constant=', &
2295                                    km_constant, ' < 0.0'
2296          CALL local_stop
2297       ELSE
2298          IF ( prandtl_number < 0.0 )  THEN
2299             IF ( myid == 0 )  PRINT*,'+++ check_parameters:  prandtl_number=',&
2300                                      prandtl_number, ' < 0.0'
2301             CALL local_stop
2302          ENDIF
2303          constant_diffusion = .TRUE.
2304
2305          IF ( prandtl_layer )  THEN
2306             IF ( myid == 0 )  PRINT*, '+++ check_parameters:  prandtl_layer ',&
2307                                       'is not allowed with fixed value of km'
2308             CALL local_stop
2309          ENDIF
2310       ENDIF
2311    ENDIF
2312
2313!
2314!-- In case of non-cyclic lateral boundaries, set the default maximum value
2315!-- for the horizontal diffusivity used within the outflow damping layer,
2316!-- and check/set the width of the damping layer
2317    IF ( bc_lr /= 'cyclic' ) THEN
2318       IF ( km_damp_max == -1.0 )  THEN
2319          km_damp_max = 0.5 * dx
2320       ENDIF
2321       IF ( outflow_damping_width == -1.0 )  THEN
2322          outflow_damping_width = MIN( 20, nx/2 )
2323       ENDIF
2324       IF ( outflow_damping_width <= 0  .OR.  outflow_damping_width > nx )  THEN
2325          IF ( myid == 0 )  PRINT*, '+++ check_parameters: outflow_damping w',&
2326                                    'idth out of range'
2327          CALL local_stop
2328       ENDIF
2329    ENDIF
2330
2331    IF ( bc_ns /= 'cyclic' )  THEN
2332       IF ( km_damp_max == -1.0 )  THEN
2333          km_damp_max = 0.5 * dy
2334       ENDIF
2335       IF ( outflow_damping_width == -1.0 )  THEN
2336          outflow_damping_width = MIN( 20, ny/2 )
2337       ENDIF
2338       IF ( outflow_damping_width <= 0  .OR.  outflow_damping_width > ny )  THEN
2339          IF ( myid == 0 )  PRINT*, '+++ check_parameters: outflow_damping w',&
2340                                    'idth out of range'
2341          CALL local_stop
2342       ENDIF
2343    ENDIF
2344
2345!
2346!-- Check value range for rif
2347    IF ( rif_min >= rif_max )  THEN
2348       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  rif_min=', rif_min, &
2349                                 ' must be less than rif_max=', rif_max
2350       CALL local_stop
2351    ENDIF
2352
2353!
2354!-- Determine upper and lower hight level indices for random perturbations
2355    IF ( disturbance_level_b == -1.0 )  THEN
2356       disturbance_level_b     = zu(nzb+3)
2357       disturbance_level_ind_b = nzb + 3
2358    ELSEIF ( disturbance_level_b < zu(3) )  THEN
2359       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_b=',&
2360                                 disturbance_level_b, ' must be >= ',zu(3),    &
2361                                 '(zu(3))'
2362       CALL local_stop
2363    ELSEIF ( disturbance_level_b > zu(nzt-2) )  THEN
2364       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_b=',&
2365                                 disturbance_level_b, ' must be <= ',zu(nzt-2),&
2366                                 '(zu(nzt-2))'
2367       CALL local_stop
2368    ELSE
2369       DO  k = 3, nzt-2
2370          IF ( disturbance_level_b <= zu(k) )  THEN
2371             disturbance_level_ind_b = k
2372             EXIT
2373          ENDIF
2374       ENDDO
2375    ENDIF
2376
2377    IF ( disturbance_level_t == -1.0 )  THEN
2378       disturbance_level_t     = zu(nzt/3)
2379       disturbance_level_ind_t = nzt / 3
2380    ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
2381       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_t=',&
2382                                 disturbance_level_t, ' must be <= ',zu(nzt-2),&
2383                                 '(zu(nzt-2))'
2384       CALL local_stop
2385    ELSEIF ( disturbance_level_t < disturbance_level_b )  THEN
2386       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_t=',&
2387                                 disturbance_level_t, ' must be >= ',          &
2388                                 'disturbance_level_b=', disturbance_level_b
2389       CALL local_stop
2390    ELSE
2391       DO  k = 3, nzt-2
2392          IF ( disturbance_level_t <= zu(k) )  THEN
2393             disturbance_level_ind_t = k
2394             EXIT
2395          ENDIF
2396       ENDDO
2397    ENDIF
2398
2399!
2400!-- Check again whether the levels determined this way are ok.
2401!-- Error may occur at automatic determination and too few grid points in
2402!-- z-direction.
2403    IF ( disturbance_level_ind_t < disturbance_level_ind_b )  THEN
2404       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  ',                &
2405                                 'disturbance_level_ind_t=',               &
2406                                 disturbance_level_ind_t, ' must be >= ',  &
2407                                 'disturbance_level_ind_b=',               &
2408                                 disturbance_level_ind_b
2409       CALL local_stop
2410    ENDIF
2411
2412!
2413!-- Determine the horizontal index range for random perturbations.
2414!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed
2415!-- near the inflow and the perturbation area is further limited to ...(1)
2416!-- after the initial phase of the flow.
2417    dist_nxl = 0;  dist_nxr = nx
2418    dist_nys = 0;  dist_nyn = ny
2419    IF ( bc_lr /= 'cyclic' )  THEN
2420       IF ( inflow_disturbance_begin == -1 )  THEN
2421          inflow_disturbance_begin = MIN( 10, nx/2 )
2422       ENDIF
2423       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )&
2424       THEN
2425          IF ( myid == 0 )  PRINT*, '+++ check_parameters: inflow_disturbance',&
2426                                    '_begin out of range'
2427          CALL local_stop
2428       ENDIF
2429       IF ( inflow_disturbance_end == -1 )  THEN
2430          inflow_disturbance_end = MIN( 100, 3*nx/4 )
2431       ENDIF
2432       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )    &
2433       THEN
2434          IF ( myid == 0 )  PRINT*, '+++ check_parameters: inflow_disturbance',&
2435                                    '_end out of range'
2436          CALL local_stop
2437       ENDIF
2438    ELSEIF ( bc_ns /= 'cyclic' )  THEN
2439       IF ( inflow_disturbance_begin == -1 )  THEN
2440          inflow_disturbance_begin = MIN( 10, ny/2 )
2441       ENDIF
2442       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )&
2443       THEN
2444          IF ( myid == 0 )  PRINT*, '+++ check_parameters: inflow_disturbance',&
2445                                    '_begin out of range'
2446          CALL local_stop
2447       ENDIF
2448       IF ( inflow_disturbance_end == -1 )  THEN
2449          inflow_disturbance_end = MIN( 100, 3*ny/4 )
2450       ENDIF
2451       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )    &
2452       THEN
2453          IF ( myid == 0 )  PRINT*, '+++ check_parameters: inflow_disturbance',&
2454                                    '_end out of range'
2455          CALL local_stop
2456       ENDIF
2457    ENDIF
2458
2459    IF ( bc_lr == 'radiation/dirichlet' )  THEN
2460       dist_nxr    = nx - inflow_disturbance_begin
2461       dist_nxl(1) = nx - inflow_disturbance_end
2462    ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
2463       dist_nxl    = inflow_disturbance_begin
2464       dist_nxr(1) = inflow_disturbance_end
2465    ENDIF
2466    IF ( bc_ns == 'dirichlet/radiation' )  THEN
2467       dist_nyn    = ny - inflow_disturbance_begin
2468       dist_nys(1) = ny - inflow_disturbance_end
2469    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
2470       dist_nys    = inflow_disturbance_begin
2471       dist_nyn(1) = inflow_disturbance_end
2472    ENDIF
2473
2474!
2475!-- Check random generator
2476    IF ( random_generator /= 'system-specific'  .AND. &
2477         random_generator /= 'numerical-recipes' )  THEN
2478       IF ( myid == 0 )  THEN
2479          PRINT*, '+++ check_parameters:'
2480          PRINT*, '    unknown random generator: random_generator=', &
2481                  random_generator
2482       ENDIF
2483       CALL local_stop
2484    ENDIF
2485
2486!
2487!-- Determine damping level index for 1D model
2488    IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
2489       IF ( damp_level_1d == -1.0 )  THEN
2490          damp_level_1d     = zu(nzt+1)
2491          damp_level_ind_1d = nzt + 1
2492       ELSEIF ( damp_level_1d < 0.0  .OR.  damp_level_1d > zu(nzt+1) )  THEN
2493          IF ( myid == 0 )  PRINT*, '+++ check_parameters:  damp_level_1d=', &
2494                                    damp_level_1d, ' must be > 0.0 and < ',  &
2495                                    'zu(nzt+1)', zu(nzt+1)
2496          CALL local_stop
2497       ELSE
2498          DO  k = 1, nzt+1
2499             IF ( damp_level_1d <= zu(k) )  THEN
2500                damp_level_ind_1d = k
2501                EXIT
2502             ENDIF
2503          ENDDO
2504       ENDIF
2505    ENDIF
2506!
2507!-- Check some other 1d-model parameters
2508    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND. &
2509         TRIM( mixing_length_1d ) /= 'blackadar' )  THEN
2510       IF ( myid == 0 )  PRINT*, '+++ check_parameters: mixing_length_1d = "', &
2511                                 TRIM( mixing_length_1d ), '" is unknown'
2512       CALL local_stop
2513    ENDIF
2514    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND. &
2515         TRIM( dissipation_1d ) /= 'detering' )  THEN
2516       IF ( myid == 0 )  PRINT*, '+++ check_parameters: dissipation_1d = "', &
2517                                 TRIM( dissipation_1d ), '" is unknown'
2518       CALL local_stop
2519    ENDIF
2520
2521!
2522!-- Set time for the next user defined restart (time_restart is the
2523!-- internal parameter for steering restart events)
2524    IF ( restart_time /= 9999999.9 )  THEN
2525       IF ( restart_time > simulated_time )  time_restart = restart_time
2526    ELSE
2527!
2528!--    In case of a restart run, set internal parameter to default (no restart)
2529!--    if the NAMELIST-parameter restart_time is omitted
2530       time_restart = 9999999.9
2531    ENDIF
2532
2533!
2534!-- Set default value of the time needed to terminate a model run
2535    IF ( termination_time_needed == -1.0 )  THEN
2536       IF ( host(1:3) == 'ibm' )  THEN
2537          termination_time_needed = 300.0
2538       ELSE
2539          termination_time_needed = 35.0
2540       ENDIF
2541    ENDIF
2542
2543!
2544!-- Check the time needed to terminate a model run
2545    IF ( host(1:3) == 't3e' )  THEN
2546!
2547!--    Time needed must be at least 30 seconds on all CRAY machines, because
2548!--    MPP_TREMAIN gives the remaining CPU time only in steps of 30 seconds
2549       IF ( termination_time_needed <= 30.0 )  THEN
2550          IF ( myid == 0 )  THEN
2551             PRINT*, '+++ check_parameters:  termination_time_needed', &
2552                      termination_time_needed
2553             PRINT*, '                       must be > 30.0 on host "', host, &
2554                     '"'
2555          ENDIF
2556          CALL local_stop
2557       ENDIF
2558    ELSEIF ( host(1:3) == 'ibm' )  THEN
2559!
2560!--    On IBM-regatta machines the time should be at least 300 seconds,
2561!--    because the job time consumed before executing palm (for compiling,
2562!--    copying of files, etc.) has to be regarded
2563       IF ( termination_time_needed < 300.0 )  THEN
2564          IF ( myid == 0 )  THEN
2565             PRINT*, '+++ WARNING: check_parameters:  termination_time_',  &
2566                         'needed', termination_time_needed
2567             PRINT*, '                                should be >= 300.0', &
2568                         ' on host "', host, '"'
2569          ENDIF
2570       ENDIF
2571    ENDIF
2572
2573
2574 END SUBROUTINE check_parameters
Note: See TracBrowser for help on using the repository browser.