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

Last change on this file since 688 was 688, checked in by gryschka, 13 years ago

Bugifx for some logical expressions (syntax was not compatible with all compilers)

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