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

Last change on this file since 150 was 147, checked in by raasch, 16 years ago

further updates for turbulent inflow: reading input data of a precursor run using a smaller total domain is working

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