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

Last change on this file since 667 was 667, checked in by suehring, 11 years ago

summary:


Gryschka:

  • Coupling with different resolution and different numbers of PEs in ocean and atmosphere is available
  • Exchange of u and v from ocean surface to atmosphere surface
  • Mirror boundary condition for u and v at the bottom are replaced by dirichlet boundary conditions
  • Inflow turbulence is now defined by flucuations around spanwise mean
  • Bugfixes for cyclic_fill and constant_volume_flow

Suehring:

  • New advection added ( Wicker and Skamarock 5th order ), therefore:
    • New module advec_ws.f90
    • Modified exchange of ghost boundaries.
    • Modified evaluation of turbulent fluxes
    • New index bounds nxlg, nxrg, nysg, nyng

advec_ws.f90


Advection scheme for scalars and momentum using the flux formulation of
Wicker and Skamarock 5th order.
Additionally the module contains of a routine using for initialisation and
steering of the statical evaluation. The computation of turbulent fluxes takes
place inside the advection routines.
In case of vector architectures Dirichlet and Radiation boundary conditions are
outstanding and not available. Furthermore simulations within topography are
not possible so far. A further routine local_diss_ij is available and is used
if a control of dissipative fluxes is desired.

check_parameters.f90


Exchange of parameters between ocean and atmosphere via PE0
Check for illegal combination of ws-scheme and timestep scheme.
Check for topography and ws-scheme.
Check for not cyclic boundary conditions in combination with ws-scheme and
loop_optimization = 'vector'.
Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.

Different processor/grid topology in atmosphere and ocean is now allowed!
Bugfixes in checking for conserve_volume_flow_mode.

exchange_horiz.f90


Dynamic exchange of ghost points with nbgp_local to ensure that no useless
ghost points exchanged in case of multigrid. type_yz(0) and type_xz(0) used for
normal grid, the remaining types used for the several grid levels.
Exchange is done via MPI-Vectors with a dynamic value of ghost points which
depend on the advection scheme. Exchange of left and right PEs is 10% faster
with MPI-Vectors than without.

flow_statistics.f90


When advection is computed with ws-scheme, turbulent fluxes are already
computed in the respective advection routines and buffered in arrays
sums_xxxx_ws_l(). This is due to a consistent treatment of statistics
with the numerics and to avoid unphysical kinks near the surface. So some if-
requests has to be done to dicern between fluxes from ws-scheme other advection
schemes. Furthermore the computation of z_i is only done if the heat flux
exceeds a minimum value. This affects only simulations of a neutral boundary
layer and is due to reasons of computations in the advection scheme.

inflow_turbulence.f90


Using nbgp recycling planes for a better resolution of the turbulent flow near
the inflow.

init_grid.f90


Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
Furthermore the allocation of arrays and steering of loops is done with these
parameters. Call of exchange_horiz are modified.
In case of dirichlet bounday condition at the bottom zu(0)=0.0
dzu_mg has to be set explicitly for a equally spaced grid near bottom.
ddzu_pres added to use a equally spaced grid near bottom.

init_pegrid.f90


Moved determination of target_id's from init_coupling
Determination of parameters needed for coupling (coupling_topology, ngp_a, ngp_o)
with different grid/processor-topology in ocean and atmosphere

Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
This distinction is due to reasons of data exchange and performance for the
normal grid and grids in poismg.
The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
New MPI-Vectors for data exchange between left and right boundaries added.
This is due to reasons of performance (10% faster).

ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values

parin.f90


Steering parameter dissipation_control added in inipar.

Makefile


Module advec_ws added.

Modules


Removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc

For coupling with different resolution in ocean and atmophere:
+nx_a, +nx_o, ny_a, +ny_o, ngp_a, ngp_o, +total_2d_o, +total_2d_a,
+coupling_topology

Buffer arrays for the left sided advective fluxes added in arrays_3d.
+flux_s_u, +flux_s_v, +flux_s_w, +diss_s_u, +diss_s_v, +diss_s_w,
+flux_s_pt, +diss_s_pt, +flux_s_e, +diss_s_e, +flux_s_q, +diss_s_q,
+flux_s_sa, +diss_s_sa
3d arrays for dissipation control added. (only necessary for vector arch.)
+var_x, +var_y, +var_z, +gamma_x, +gamma_y, +gamma_z
Default of momentum_advec and scalar_advec changed to 'ws-scheme' .
+exchange_mg added in control_parameters to steer the data exchange.
Parameters +nbgp, +nxlg, +nxrg, +nysg, +nyng added in indices.
flag array +boundary_flags added in indices to steer the degradation of order
of the advective fluxes when non-cyclic boundaries are used.
MPI-datatypes +type_y, +type_y_int and +type_yz for data_exchange added in
pegrid.
+sums_wsus_ws_l, +sums_wsvs_ws_l, +sums_us2_ws_l, +sums_vs2_ws_l,
+sums_ws2_ws_l, +sums_wspts_ws_l, +sums_wssas_ws_l, +sums_wsqs_ws_l
and +weight_substep added in statistics to steer the statistical evaluation
of turbulent fluxes in the advection routines.
LOGICALS +ws_scheme_sca and +ws_scheme_mom added to get a better performance
in prognostic_equations.
LOGICAL +dissipation_control control added to steer numerical dissipation
in ws-scheme.

Changed length of string run_description_header

pres.f90


New allocation of tend when ws-scheme and multigrid is used. This is due to
reasons of perforance of the data_exchange. The same is done with p after
poismg is called.
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
multigrid is used. Calls of exchange_horiz are modified.

bugfix: After pressure correction no volume flow correction in case of
non-cyclic boundary conditions
(has to be done only before pressure correction)

Call of SOR routine is referenced with ddzu_pres.

prognostic_equations.f90


Calls of the advection routines with WS5 added.
Calls of ws_statistics added to set the statistical arrays to zero after each
time step.

advec_particles.f90


Declaration of de_dx, de_dy, de_dz adapted to additional ghost points.
Furthermore the calls of exchange_horiz were modified.

asselin_filter.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

average_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

boundary_conds.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
Removed mirror boundary conditions for u and v at the bottom in case of
ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
in init_3d_model

calc_liquid_water_content.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

calc_spectra.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for
allocation of tend.

check_open.f90


Output of total array size was adapted to nbgp.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays local_2d and total_2d.
Calls of exchange_horiz are modified.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Skip-value skip_do_avs changed to a dynamic adaption of ghost points.

data_output_mask.f90


Calls of exchange_horiz are modified.

diffusion_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_s.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_u.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_v.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_w.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

exchange_horiz_2d.f90


Dynamic exchange of ghost points with nbgp, which depends on the advection
scheme. Exchange between left and right PEs is now done with MPI-vectors.

global_min_max.f90


Adapting of the index arrays, because MINLOC assumes lowerbound
at 1 and not at nbgp.

init_3d_model.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Call ws_init to initialize arrays needed for statistical evaluation and
optimization when ws-scheme is used.
Initial volume flow is now calculated by using the variable hom_sum.
Therefore the correction of initial volume flow for non-flat topography
removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
representative for the height z0

Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)

init_coupling.f90


determination of target_id's moved to init_pegrid

init_pt_anomaly.f90


Call of exchange_horiz are modified.

init_rankine.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

init_slope.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

header.f90


Output of advection scheme.

poismg.f90


Calls of exchange_horiz are modified.

prandtl_fluxes.f90


Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

production_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

read_3d_binary.f90


+/- 1 replaced with +/- nbgp when swapping and allocating variables.

sor.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Call of exchange_horiz are modified.
bug removed in declaration of ddzw(), nz replaced by nzt+1

subsidence.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

sum_up_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

surface_coupler.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in
MPI_SEND() and MPI_RECV.
additional case for nonequivalent processor and grid topopolgy in ocean and
atmosphere added (coupling_topology = 1)

Added exchange of u and v from Ocean to Atmosphere

time_integration.f90


Calls of exchange_horiz are modified.
Adaption to slooping surface.

timestep.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_3d_data_averaging.f90, user_data_output_2d.f90, user_data_output_3d.f90,
user_actions.f90, user_init.f90, user_init_plant_canopy.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_read_restart_data.f90


Allocation with nbgp.

wall_fluxes.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

write_compressed.f90


Array bounds and nx, ny adapted with nbgp.

sor.f90


bug removed in declaration of ddzw(), nz replaced by nzt+1

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