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

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

WS5 is available in combination with topography. Version number changed from 3.8 to 3.8a. Modification in init_pt_anomaly.

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