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

Last change on this file since 212 was 206, checked in by raasch, 15 years ago

ocean-atmosphere coupling realized with MPI-1, adjustments in mrun, mbuild, subjob for lcxt4

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