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

Last change on this file since 177 was 177, checked in by steinfeld, 16 years ago

Bug in the construction of initial vertical profiles fixed

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