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

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

last commit documented

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