source: palm/tags/release-3.4a/SOURCE/check_parameters.f90 @ 247

Last change on this file since 247 was 139, checked in by raasch, 16 years ago

New:
---

Plant canopy model of Watanabe (2004,BLM 112,307-341) added.
It can be switched on by the inipar parameter plant_canopy.
The inipar parameter canopy_mode can be used to prescribe a
plant canopy type. The default case is a homogeneous plant
canopy. Heterogeneous distributions of the leaf area
density and the canopy drag coefficient can be defined in the
new routine user_init_plant_canopy (user_interface).
The inipar parameters lad_surface, lad_vertical_gradient and
lad_vertical_gradient_level can be used in order to
prescribe the vertical profile of leaf area density. The
inipar parameter drag_coefficient determines the canopy
drag coefficient.
Finally, the inipar parameter pch_index determines the
index of the upper boundary of the plant canopy.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

Pressure boundary conditions for vertical walls added to the multigrid solver.
They are applied using new wall flag arrays (wall_flags_..) which are defined
for each grid level. New argument gls added to routine user_init_grid
(user_interface).

Frequence of sorting particles can be controlled with new particles_par
parameter dt_sort_particles. Sorting is moved from the SGS timestep loop in
advec_particles after the end of this loop.

advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, parin, plant_canopy_model, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary

Changed:


Redefine initial nzb_local as the actual total size of topography (later the
extent of topography in nzb_local is reduced by 1dx at the E topography walls
and by 1dy at the N topography walls to form the basis for nzb_s_inner);
for consistency redefine 'single_building' case.

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered velocity
components and their products, procucts of scalars and velocity components),
respectively.

Allow two instead of one digit to specify isosurface and slicer variables.

Status of 3D-volume NetCDF data file only depends on switch netcdf_64bit_3d (check_open)

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
and wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

  • it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

  • it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

buoyancy, check_open, data_output_dvrp, diffusion_s, diffusivities, flow_statistics, header, init_3d_model, init_dvrp, init_grid, modules, prognostic_equations

Errors:


Bugfix: summation of sums_l_l in diffusivities.

Several bugfixes in the ocean part: Initial density rho is calculated
(init_ocean). Error in initializing u_init and v_init removed
(check_parameters). Calculation of density flux now starts from
nzb+1 (production_e).

Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
numbers along y, small bugfixes in the SGS part (advec_particles)

Bugfix: model_string needed a default value (combine_plot_fields)

Bugfix: wavenumber calculation for even nx in routines maketri (poisfft)

Bugfix: assignment of fluxes at walls

Bugfix: absolute value of f must be used when calculating the Blackadar mixing length (init_1d_model)

advec_particles, check_parameters, combine_plot_fields, diffusion_s, diffusivities, init_ocean, init_1d_model, poisfft, production_e

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