source: palm/tags/release-3.7/SOURCE/check_parameters.f90 @ 3806

Last change on this file since 3806 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:


lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:


Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

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