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

Last change on this file since 760 was 708, checked in by raasch, 13 years ago

last commit documented

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