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

Last change on this file since 1683 was 1683, checked in by knoop, 9 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 175.0 KB
Line 
1!> @file check_parameters.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2014 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: check_parameters.f90 1683 2015-10-07 23:57:51Z knoop $
26!
27! 1682 2015-10-07 23:56:08Z knoop
28! Code annotations made doxygen readable
29!
30! 1606 2015-06-29 10:43:37Z maronga
31! Added check for use of RRTMG without netCDF.
32!
33! 1587 2015-05-04 14:19:01Z maronga
34! Added check for set of albedos when using RRTMG
35!
36! 1585 2015-04-30 07:05:52Z maronga
37! Added support for RRTMG
38!
39! 1575 2015-03-27 09:56:27Z raasch
40! multigrid_fast added as allowed pressure solver
41!
42! 1573 2015-03-23 17:53:37Z suehring
43! Bugfix: check for advection schemes in case of non-cyclic boundary conditions
44!
45! 1557 2015-03-05 16:43:04Z suehring
46! Added checks for monotonic limiter
47!
48! 1555 2015-03-04 17:44:27Z maronga
49! Added output of r_a and r_s. Renumbering of LSM PA-messages.
50!
51! 1553 2015-03-03 17:33:54Z maronga
52! Removed check for missing soil temperature profile as default values were added.
53!
54! 1551 2015-03-03 14:18:16Z maronga
55! Added various checks for land surface and radiation model. In the course of this
56! action, the length of the variable var has to be increased
57!
58! 1504 2014-12-04 09:23:49Z maronga
59! Bugfix: check for large_scale forcing got lost.
60!
61! 1500 2014-12-03 17:42:41Z maronga
62! Boundary conditions changed to dirichlet for land surface model
63!
64! 1496 2014-12-02 17:25:50Z maronga
65! Added checks for the land surface model
66!
67! 1484 2014-10-21 10:53:05Z kanani
68! Changes due to new module structure of the plant canopy model:
69!   module plant_canopy_model_mod added,
70!   checks regarding usage of new method for leaf area density profile
71!   construction added,
72!   lad-profile construction moved to new subroutine init_plant_canopy within
73!   the module plant_canopy_model_mod,
74!   drag_coefficient renamed to canopy_drag_coeff.
75! Missing KIND-attribute for REAL constant added
76!
77! 1455 2014-08-29 10:47:47Z heinze
78! empty time records in volume, cross-section and masked data output prevented 
79! in case of non-parallel netcdf-output in restart runs
80!
81! 1429 2014-07-15 12:53:45Z knoop
82! run_description_header exended to provide ensemble_member_nr if specified
83!
84! 1425 2014-07-05 10:57:53Z knoop
85! bugfix: perturbation domain modified for parallel random number generator
86!
87! 1402 2014-05-09 14:25:13Z raasch
88! location messages modified
89!
90! 1400 2014-05-09 14:03:54Z knoop
91! Check random generator extended by option random-parallel
92!
93! 1384 2014-05-02 14:31:06Z raasch
94! location messages added
95!
96! 1365 2014-04-22 15:03:56Z boeske
97! Usage of large scale forcing for pt and q enabled:
98! output for profiles of large scale advection (td_lsa_lpt, td_lsa_q),
99! large scale subsidence (td_sub_lpt, td_sub_q)
100! and nudging tendencies (td_nud_lpt, td_nud_q, td_nud_u and td_nud_v) added,
101! check if use_subsidence_tendencies is used correctly
102!
103! 1361 2014-04-16 15:17:48Z hoffmann
104! PA0363 removed
105! PA0362 changed
106!
107! 1359 2014-04-11 17:15:14Z hoffmann
108! Do not allow the execution of PALM with use_particle_tails, since particle
109! tails are currently not supported by our new particle structure.
110!
111! PA0084 not necessary for new particle structure
112!
113! 1353 2014-04-08 15:21:23Z heinze
114! REAL constants provided with KIND-attribute
115!
116! 1330 2014-03-24 17:29:32Z suehring
117! In case of SGS-particle velocity advection of TKE is also allowed with
118! dissipative 5th-order scheme.
119!
120! 1327 2014-03-21 11:00:16Z raasch
121! "baroclinicity" renamed "baroclinity", "ocean version" replaced by "ocean mode"
122! bugfix: duplicate error message 56 removed,
123! check of data_output_format and do3d_compress removed
124!
125! 1322 2014-03-20 16:38:49Z raasch
126! some REAL constants defined as wp-kind
127!
128! 1320 2014-03-20 08:40:49Z raasch
129! Kind-parameters added to all INTEGER and REAL declaration statements,
130! kinds are defined in new module kinds,
131! revision history before 2012 removed,
132! comment fields (!:) to be used for variable explanations added to
133! all variable declaration statements
134!
135! 1308 2014-03-13 14:58:42Z fricke
136! +netcdf_data_format_save
137! Calculate fixed number of output time levels for parallel netcdf output.
138! For masked data, parallel netcdf output is not tested so far, hence
139! netcdf_data_format is switched back to non-paralell output.
140!
141! 1299 2014-03-06 13:15:21Z heinze
142! enable usage of large_scale subsidence in combination with large_scale_forcing
143! output for profile of large scale vertical velocity w_subs added
144!
145! 1276 2014-01-15 13:40:41Z heinze
146! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
147!
148! 1241 2013-10-30 11:36:58Z heinze
149! output for profiles of ug and vg added
150! set surface_heatflux and surface_waterflux also in dependence on
151! large_scale_forcing
152! checks for nudging and large scale forcing from external file
153!
154! 1236 2013-09-27 07:21:13Z raasch
155! check number of spectra levels
156!
157! 1216 2013-08-26 09:31:42Z raasch
158! check for transpose_compute_overlap (temporary)
159!
160! 1214 2013-08-21 12:29:17Z kanani
161! additional check for simultaneous use of vertical grid stretching
162! and particle advection
163!
164! 1212 2013-08-15 08:46:27Z raasch
165! checks for poisfft_hybrid removed
166!
167! 1210 2013-08-14 10:58:20Z raasch
168! check for fftw
169!
170! 1179 2013-06-14 05:57:58Z raasch
171! checks and settings of buoyancy parameters and switches revised,
172! initial profile for rho added to hom (id=77)
173!
174! 1174 2013-05-31 10:28:08Z gryschka
175! Bugfix in computing initial profiles for ug, vg, lad, q in case of Atmosphere
176!
177! 1159 2013-05-21 11:58:22Z fricke
178! bc_lr/ns_dirneu/neudir removed
179!
180! 1115 2013-03-26 18:16:16Z hoffmann
181! unused variables removed
182! drizzle can be used without precipitation
183!
184! 1111 2013-03-08 23:54:10Z raasch
185! ibc_p_b = 2 removed
186!
187! 1103 2013-02-20 02:15:53Z raasch
188! Bugfix: turbulent inflow must not require cyclic fill in restart runs
189!
190! 1092 2013-02-02 11:24:22Z raasch
191! unused variables removed
192!
193! 1069 2012-11-28 16:18:43Z maronga
194! allow usage of topography in combination with cloud physics
195!
196! 1065 2012-11-22 17:42:36Z hoffmann
197! Bugfix: It is not allowed to use cloud_scheme = seifert_beheng without
198!         precipitation in order to save computational resources.
199!
200! 1060 2012-11-21 07:19:51Z raasch
201! additional check for parameter turbulent_inflow
202!
203! 1053 2012-11-13 17:11:03Z hoffmann
204! necessary changes for the new two-moment cloud physics scheme added:
205! - check cloud physics scheme (Kessler or Seifert and Beheng)
206! - plant_canopy is not allowed
207! - currently, only cache loop_optimization is allowed
208! - initial profiles of nr, qr
209! - boundary condition of nr, qr
210! - check output quantities (qr, nr, prr)
211!
212! 1036 2012-10-22 13:43:42Z raasch
213! code put under GPL (PALM 3.9)
214!
215! 1031/1034 2012-10-22 11:32:49Z raasch
216! check of netcdf4 parallel file support
217!
218! 1019 2012-09-28 06:46:45Z raasch
219! non-optimized version of prognostic_equations not allowed any more
220!
221! 1015 2012-09-27 09:23:24Z raasch
222! acc allowed for loop optimization,
223! checks for adjustment of mixing length to the Prandtl mixing length removed
224!
225! 1003 2012-09-14 14:35:53Z raasch
226! checks for cases with unequal subdomain sizes removed
227!
228! 1001 2012-09-13 14:08:46Z raasch
229! all actions concerning leapfrog- and upstream-spline-scheme removed
230!
231! 996 2012-09-07 10:41:47Z raasch
232! little reformatting
233
234! 978 2012-08-09 08:28:32Z fricke
235! setting of bc_lr/ns_dirneu/neudir
236! outflow damping layer removed
237! check for z0h*
238! check for pt_damping_width
239!
240! 964 2012-07-26 09:14:24Z raasch
241! check of old profil-parameters removed
242!
243! 940 2012-07-09 14:31:00Z raasch
244! checks for parameter neutral
245!
246! 924 2012-06-06 07:44:41Z maronga
247! Bugfix: preprocessor directives caused error during compilation
248!
249! 892 2012-05-02 13:51:44Z maronga
250! Bugfix for parameter file check ( excluding __netcdf4 )
251!
252! 866 2012-03-28 06:44:41Z raasch
253! use only 60% of the geostrophic wind as translation speed in case of Galilean
254! transformation and use_ug_for_galilei_tr = .T. in order to mimimize the
255! timestep
256!
257! 861 2012-03-26 14:18:34Z suehring
258! Check for topography and ws-scheme removed.
259! Check for loop_optimization = 'vector' and ws-scheme removed.
260!
261! 845 2012-03-07 10:23:05Z maronga
262! Bugfix: exclude __netcdf4 directive part from namelist file check compilation
263!
264! 828 2012-02-21 12:00:36Z raasch
265! check of collision_kernel extended
266!
267! 825 2012-02-19 03:03:44Z raasch
268! check for collision_kernel and curvature_solution_effects
269!
270! 809 2012-01-30 13:32:58Z maronga
271! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
272!
273! 807 2012-01-25 11:53:51Z maronga
274! New cpp directive "__check" implemented which is used by check_namelist_files
275!
276! Revision 1.1  1997/08/26 06:29:23  raasch
277! Initial revision
278!
279!
280! Description:
281! ------------
282!> Check control parameters and deduce further quantities.
283!------------------------------------------------------------------------------!
284 SUBROUTINE check_parameters
285 
286
287    USE arrays_3d
288    USE cloud_parameters
289    USE constants
290    USE control_parameters
291    USE dvrp_variables
292    USE grid_variables
293    USE indices
294    USE land_surface_model_mod
295    USE kinds
296    USE model_1d
297    USE netcdf_control
298    USE particle_attributes
299    USE pegrid
300    USE plant_canopy_model_mod
301    USE profil_parameter
302    USE radiation_model_mod
303    USE spectrum
304    USE statistics
305    USE subsidence_mod
306    USE statistics
307    USE transpose_indices
308
309    IMPLICIT NONE
310
311    CHARACTER (LEN=1)   ::  sq                       !<
312    CHARACTER (LEN=15)  ::  var                      !<
313    CHARACTER (LEN=7)   ::  unit                     !<
314    CHARACTER (LEN=8)   ::  date                     !<
315    CHARACTER (LEN=10)  ::  time                     !<
316    CHARACTER (LEN=40)  ::  coupling_string          !<
317    CHARACTER (LEN=100) ::  action                   !<
318
319    INTEGER(iwp) ::  i                               !<
320    INTEGER(iwp) ::  ilen                            !<
321    INTEGER(iwp) ::  iremote = 0                     !<
322    INTEGER(iwp) ::  j                               !<
323    INTEGER(iwp) ::  k                               !<
324    INTEGER(iwp) ::  kk                              !<
325    INTEGER(iwp) ::  netcdf_data_format_save         !<
326    INTEGER(iwp) ::  position                        !<
327    INTEGER(iwp) ::  prec                            !<
328   
329    LOGICAL     ::  found                            !<
330    LOGICAL     ::  ldum                             !<
331   
332    REAL(wp)    ::  gradient                         !<
333    REAL(wp)    ::  remote = 0.0_wp                  !<
334    REAL(wp)    ::  simulation_time_since_reference  !<
335
336
337    CALL location_message( 'checking parameters', .FALSE. )
338
339!
340!-- Check for overlap combinations, which are not realized yet
341    IF ( transpose_compute_overlap )  THEN
342       IF ( numprocs == 1 )  STOP '+++ transpose-compute-overlap not implemented for single PE runs'
343#if defined( __openacc )
344       STOP '+++ transpose-compute-overlap not implemented for GPU usage'
345#endif
346    ENDIF
347
348!
349!-- Warning, if host is not set
350    IF ( host(1:1) == ' ' )  THEN
351       message_string = '"host" is not set. Please check that environment ' // &
352                        'variable "localhost" & is set before running PALM'
353       CALL message( 'check_parameters', 'PA0001', 0, 0, 0, 6, 0 )
354    ENDIF
355
356!
357!-- Check the coupling mode
358    IF ( coupling_mode /= 'uncoupled'            .AND.  &
359         coupling_mode /= 'atmosphere_to_ocean'  .AND.  &
360         coupling_mode /= 'ocean_to_atmosphere' )  THEN
361       message_string = 'illegal coupling mode: ' // TRIM( coupling_mode )
362       CALL message( 'check_parameters', 'PA0002', 1, 2, 0, 6, 0 )
363    ENDIF
364
365!
366!-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
367    IF ( coupling_mode /= 'uncoupled')  THEN
368
369       IF ( dt_coupling == 9999999.9_wp )  THEN
370          message_string = 'dt_coupling is not set but required for coup' // &
371                           'ling mode "' //  TRIM( coupling_mode ) // '"'
372          CALL message( 'check_parameters', 'PA0003', 1, 2, 0, 6, 0 )
373       ENDIF
374
375#if defined( __parallel )
376
377!
378!--    NOTE: coupled runs have not been implemented in the check_namelist_files
379!--    program.
380!--    check_namelist_files will need the following information of the other
381!--    model (atmosphere/ocean).
382!       dt_coupling = remote
383!       dt_max = remote
384!       restart_time = remote
385!       dt_restart= remote
386!       simulation_time_since_reference = remote
387!       dx = remote
388
389
390#if ! defined( __check )
391       IF ( myid == 0 ) THEN
392          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, &
393                         ierr )
394          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, &
395                         status, ierr )
396       ENDIF
397       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
398#endif     
399       IF ( dt_coupling /= remote )  THEN
400          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
401                 '": dt_coupling = ', dt_coupling, '& is not equal to ',       &
402                 'dt_coupling_remote = ', remote
403          CALL message( 'check_parameters', 'PA0004', 1, 2, 0, 6, 0 )
404       ENDIF
405       IF ( dt_coupling <= 0.0_wp )  THEN
406#if ! defined( __check )
407          IF ( myid == 0  ) THEN
408             CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
409             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, &
410                            status, ierr )
411          ENDIF   
412          CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
413#endif         
414          dt_coupling = MAX( dt_max, remote )
415          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
416                 '": dt_coupling <= 0.0 & is not allowed and is reset to ',    &
417                 'MAX(dt_max(A,O)) = ', dt_coupling
418          CALL message( 'check_parameters', 'PA0005', 0, 1, 0, 6, 0 )
419       ENDIF
420#if ! defined( __check )
421       IF ( myid == 0 ) THEN
422          CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
423                         ierr )
424          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, &
425                         status, ierr )
426       ENDIF
427       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
428#endif     
429       IF ( restart_time /= remote )  THEN
430          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
431                 '": restart_time = ', restart_time, '& is not equal to ',     &
432                 'restart_time_remote = ', remote
433          CALL message( 'check_parameters', 'PA0006', 1, 2, 0, 6, 0 )
434       ENDIF
435#if ! defined( __check )
436       IF ( myid == 0 ) THEN
437          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, &
438                         ierr )
439          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, &
440                         status, ierr )
441       ENDIF   
442       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
443#endif     
444       IF ( dt_restart /= remote )  THEN
445          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
446                 '": dt_restart = ', dt_restart, '& is not equal to ',         &
447                 'dt_restart_remote = ', remote
448          CALL message( 'check_parameters', 'PA0007', 1, 2, 0, 6, 0 )
449       ENDIF
450
451       simulation_time_since_reference = end_time - coupling_start_time
452#if ! defined( __check )
453       IF  ( myid == 0 ) THEN
454          CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &
455                         14, comm_inter, ierr )
456          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, &
457                         status, ierr )   
458       ENDIF
459       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
460#endif     
461       IF ( simulation_time_since_reference /= remote )  THEN
462          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
463                 '": simulation_time_since_reference = ',                      &
464                 simulation_time_since_reference, '& is not equal to ',        &
465                 'simulation_time_since_reference_remote = ', remote
466          CALL message( 'check_parameters', 'PA0008', 1, 2, 0, 6, 0 )
467       ENDIF
468
469#if ! defined( __check )
470       IF ( myid == 0 ) THEN
471          CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
472          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, &
473                                                             status, ierr )
474       ENDIF
475       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
476
477#endif
478       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
479
480          IF ( dx < remote ) THEN
481             WRITE( message_string, * ) 'coupling mode "', &
482                   TRIM( coupling_mode ),                  &
483           '": dx in Atmosphere is not equal to or not larger then dx in ocean'
484             CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
485          ENDIF
486
487          IF ( (nx_a+1)*dx /= (nx_o+1)*remote )  THEN
488             WRITE( message_string, * ) 'coupling mode "', &
489                    TRIM( coupling_mode ), &
490             '": Domain size in x-direction is not equal in ocean and atmosphere'
491             CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
492          ENDIF
493
494       ENDIF
495
496#if ! defined( __check )
497       IF ( myid == 0) THEN
498          CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
499          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, &
500                         status, ierr )
501       ENDIF
502       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
503#endif
504       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
505
506          IF ( dy < remote )  THEN
507             WRITE( message_string, * ) 'coupling mode "', &
508                    TRIM( coupling_mode ), &
509                 '": dy in Atmosphere is not equal to or not larger then dy in ocean'
510             CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
511          ENDIF
512
513          IF ( (ny_a+1)*dy /= (ny_o+1)*remote )  THEN
514             WRITE( message_string, * ) 'coupling mode "', &
515                   TRIM( coupling_mode ), &
516             '": Domain size in y-direction is not equal in ocean and atmosphere'
517             CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
518          ENDIF
519
520          IF ( MOD(nx_o+1,nx_a+1) /= 0 )  THEN
521             WRITE( message_string, * ) 'coupling mode "', &
522                   TRIM( coupling_mode ), &
523             '": nx+1 in ocean is not divisible without remainder with nx+1 in', & 
524             ' atmosphere'
525             CALL message( 'check_parameters', 'PA0339', 1, 2, 0, 6, 0 )
526          ENDIF
527
528          IF ( MOD(ny_o+1,ny_a+1) /= 0 )  THEN
529             WRITE( message_string, * ) 'coupling mode "', &
530                   TRIM( coupling_mode ), &
531             '": ny+1 in ocean is not divisible without remainder with ny+1 in', & 
532             ' atmosphere'
533             CALL message( 'check_parameters', 'PA0340', 1, 2, 0, 6, 0 )
534          ENDIF
535
536       ENDIF
537#else
538       WRITE( message_string, * ) 'coupling requires PALM to be called with', &
539            ' ''mrun -K parallel'''
540       CALL message( 'check_parameters', 'PA0141', 1, 2, 0, 6, 0 )
541#endif
542    ENDIF
543
544#if defined( __parallel ) && ! defined ( __check )
545!
546!-- Exchange via intercommunicator
547    IF ( coupling_mode == 'atmosphere_to_ocean' .AND. myid == 0 )  THEN
548       CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter, &
549                      ierr )
550    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' .AND. myid == 0)  THEN
551       CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19, &
552                      comm_inter, status, ierr )
553    ENDIF
554    CALL MPI_BCAST( humidity_remote, 1, MPI_LOGICAL, 0, comm2d, ierr)
555   
556#endif
557
558
559!
560!-- Generate the file header which is used as a header for most of PALM's
561!-- output files
562    CALL DATE_AND_TIME( date, time )
563    run_date = date(7:8)//'-'//date(5:6)//'-'//date(3:4)
564    run_time = time(1:2)//':'//time(3:4)//':'//time(5:6)
565    IF ( coupling_mode == 'uncoupled' )  THEN
566       coupling_string = ''
567    ELSEIF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
568       coupling_string = ' coupled (atmosphere)'
569    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
570       coupling_string = ' coupled (ocean)'
571    ENDIF       
572
573    IF ( ensemble_member_nr /= 0 )  THEN
574       WRITE ( run_description_header,                                         &
575                  '(A,2X,A,2X,A,A,A,I2.2,A,2X,A,I2.2,2X,A,A,2X,A,1X,A)' )      &
576              TRIM( version ), TRIM( revision ), 'run: ',                      &
577              TRIM( run_identifier ), '.', runnr, TRIM( coupling_string ),     &
578              'en-no: ', ensemble_member_nr,'host: ', TRIM( host ),            &
579              run_date, run_time
580    ELSE
581       WRITE ( run_description_header,                                         &
582                  '(A,2X,A,2X,A,A,A,I2.2,A,2X,A,A,2X,A,1X,A)' )                &
583              TRIM( version ), TRIM( revision ), 'run: ',                      &
584              TRIM( run_identifier ), '.', runnr, TRIM( coupling_string ),     &
585              'host: ', TRIM( host ), run_date, run_time
586    ENDIF
587!
588!-- Check the general loop optimization method
589    IF ( loop_optimization == 'default' )  THEN
590       IF ( host(1:3) == 'nec' )  THEN
591          loop_optimization = 'vector'
592       ELSE
593          loop_optimization = 'cache'
594       ENDIF
595    ENDIF
596
597    SELECT CASE ( TRIM( loop_optimization ) )
598
599       CASE ( 'acc', 'cache', 'vector' )
600          CONTINUE
601
602       CASE DEFAULT
603          message_string = 'illegal value given for loop_optimization: "' // &
604                           TRIM( loop_optimization ) // '"'
605          CALL message( 'check_parameters', 'PA0013', 1, 2, 0, 6, 0 )
606
607    END SELECT
608
609!
610!-- Check if vertical grid stretching is used together with particles
611    IF ( dz_stretch_level < 100000.0_wp .AND. particle_advection )  THEN
612       message_string = 'Vertical grid stretching is not allowed together ' // &
613                        'with particle advection.'
614       CALL message( 'check_parameters', 'PA0017', 1, 2, 0, 6, 0 )
615    ENDIF
616
617!
618!--
619    IF ( use_particle_tails )  THEN
620       message_string = 'Particle tails are currently not available due ' //   &
621                        'to the new particle structure.'
622       CALL message( 'check_parameters', 'PA0392', 1, 2, 0, 6, 0 )
623    ENDIF
624
625!
626!-- Check topography setting (check for illegal parameter combinations)
627    IF ( topography /= 'flat' )  THEN
628       action = ' '
629       IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme'     &
630      .AND. scalar_advec /= 'ws-scheme-mono' )  THEN
631          WRITE( action, '(A,A)' )  'scalar_advec = ', scalar_advec
632       ENDIF
633       IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' ) &
634       THEN
635          WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
636       ENDIF
637       IF ( psolver == 'sor' )  THEN
638          WRITE( action, '(A,A)' )  'psolver = ', psolver
639       ENDIF
640       IF ( sloping_surface )  THEN
641          WRITE( action, '(A)' )  'sloping surface = .TRUE.'
642       ENDIF
643       IF ( galilei_transformation )  THEN
644          WRITE( action, '(A)' )  'galilei_transformation = .TRUE.'
645       ENDIF
646       IF ( cloud_physics )  THEN
647          WRITE( action, '(A)' )  'cloud_physics = .TRUE.'
648       ENDIF
649       IF ( cloud_droplets )  THEN
650          WRITE( action, '(A)' )  'cloud_droplets = .TRUE.'
651       ENDIF
652       IF ( .NOT. prandtl_layer )  THEN
653          WRITE( action, '(A)' )  'prandtl_layer = .FALSE.'
654       ENDIF
655       IF ( action /= ' ' )  THEN
656          message_string = 'a non-flat topography does not allow ' // &
657                           TRIM( action )
658          CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
659       ENDIF
660!
661!--    In case of non-flat topography, check whether the convention how to
662!--    define the topography grid has been set correctly, or whether the default
663!--    is applicable. If this is not possible, abort.
664       IF ( TRIM( topography_grid_convention ) == ' ' )  THEN
665          IF ( TRIM( topography ) /= 'single_building' .AND.  &
666               TRIM( topography ) /= 'single_street_canyon' .AND.  &
667               TRIM( topography ) /= 'read_from_file' )  THEN
668!--          The default value is not applicable here, because it is only valid
669!--          for the two standard cases 'single_building' and 'read_from_file'
670!--          defined in init_grid.
671             WRITE( message_string, * )  &
672                  'The value for "topography_grid_convention" ',  &
673                  'is not set. Its default value is & only valid for ',  &
674                  '"topography" = ''single_building'', ',  &
675                  '''single_street_canyon'' & or ''read_from_file''.',  &
676                  ' & Choose ''cell_edge'' or ''cell_center''.'
677             CALL message( 'user_check_parameters', 'PA0239', 1, 2, 0, 6, 0 )
678          ELSE
679!--          The default value is applicable here.
680!--          Set convention according to topography.
681             IF ( TRIM( topography ) == 'single_building' .OR.  &
682                  TRIM( topography ) == 'single_street_canyon' )  THEN
683                topography_grid_convention = 'cell_edge'
684             ELSEIF ( TRIM( topography ) == 'read_from_file' )  THEN
685                topography_grid_convention = 'cell_center'
686             ENDIF
687          ENDIF
688       ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND.  &
689                TRIM( topography_grid_convention ) /= 'cell_center' )  THEN
690          WRITE( message_string, * )  &
691               'The value for "topography_grid_convention" is ', &
692               'not recognized. & Choose ''cell_edge'' or ''cell_center''.'
693          CALL message( 'user_check_parameters', 'PA0240', 1, 2, 0, 6, 0 )
694       ENDIF
695
696    ENDIF
697
698!
699!-- Check ocean setting
700    IF ( ocean )  THEN
701
702       action = ' '
703       IF ( action /= ' ' )  THEN
704          message_string = 'ocean = .T. does not allow ' // TRIM( action )
705          CALL message( 'check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
706       ENDIF
707
708    ELSEIF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.  &
709             TRIM( coupling_char ) == '_O' )  THEN
710
711!
712!--    Check whether an (uncoupled) atmospheric run has been declared as an
713!--    ocean run (this setting is done via mrun-option -y)
714
715       message_string = 'ocean = .F. does not allow coupling_char = "' // &
716                        TRIM( coupling_char ) // '" set by mrun-option "-y"'
717       CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
718
719    ENDIF
720!
721!-- Check cloud scheme
722    IF ( cloud_scheme == 'seifert_beheng' )  THEN
723       icloud_scheme = 0
724    ELSEIF ( cloud_scheme == 'kessler' )  THEN
725       icloud_scheme = 1
726    ELSE
727       message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
728                        TRIM( cloud_scheme ) // '"'
729       CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 )
730    ENDIF
731!
732!-- Check whether there are any illegal values
733!-- Pressure solver:
734    IF ( psolver /= 'poisfft'  .AND.  psolver /= 'sor'  .AND.                  &
735         psolver /= 'multigrid'  .AND.  psolver /= 'multigrid_fast' )  THEN
736       message_string = 'unknown solver for perturbation pressure: psolver' // &
737                        ' = "' // TRIM( psolver ) // '"'
738       CALL message( 'check_parameters', 'PA0016', 1, 2, 0, 6, 0 )
739    ENDIF
740
741    IF ( psolver(1:9) == 'multigrid' )  THEN
742       IF ( cycle_mg == 'w' )  THEN
743          gamma_mg = 2
744       ELSEIF ( cycle_mg == 'v' )  THEN
745          gamma_mg = 1
746       ELSE
747          message_string = 'unknown multigrid cycle: cycle_mg = "' // &
748                           TRIM( cycle_mg ) // '"'
749          CALL message( 'check_parameters', 'PA0020', 1, 2, 0, 6, 0 )
750       ENDIF
751    ENDIF
752
753    IF ( fft_method /= 'singleton-algorithm'  .AND.  &
754         fft_method /= 'temperton-algorithm'  .AND.  &
755         fft_method /= 'fftw'                 .AND.  &
756         fft_method /= 'system-specific' )  THEN
757       message_string = 'unknown fft-algorithm: fft_method = "' // &
758                        TRIM( fft_method ) // '"'
759       CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
760    ENDIF
761   
762    IF( momentum_advec == 'ws-scheme' .AND. & 
763        .NOT. call_psolver_at_all_substeps  ) THEN
764        message_string = 'psolver must be called at each RK3 substep when "'//&
765                      TRIM(momentum_advec) // ' "is used for momentum_advec'
766        CALL message( 'check_parameters', 'PA0344', 1, 2, 0, 6, 0 )
767    END IF
768!
769!-- Advection schemes:
770    IF ( momentum_advec /= 'pw-scheme'  .AND.  momentum_advec /= 'ws-scheme' ) &
771    THEN
772       message_string = 'unknown advection scheme: momentum_advec = "' // &
773                        TRIM( momentum_advec ) // '"'
774       CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
775    ENDIF
776    IF ( ( momentum_advec == 'ws-scheme' .OR.  scalar_advec == 'ws-scheme'     &
777           .OR. scalar_advec == 'ws-scheme-mono' )                             &
778           .AND. ( timestep_scheme == 'euler' .OR.                             &
779                   timestep_scheme == 'runge-kutta-2' ) )                      &
780    THEN
781       message_string = 'momentum_advec or scalar_advec = "' &
782         // TRIM( momentum_advec ) // '" is not allowed with timestep_scheme = "' // &
783         TRIM( timestep_scheme ) // '"'
784       CALL message( 'check_parameters', 'PA0023', 1, 2, 0, 6, 0 )
785    ENDIF
786    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' .AND. &
787         scalar_advec /= 'ws-scheme-mono' .AND. scalar_advec /= 'bc-scheme' )  &
788    THEN
789       message_string = 'unknown advection scheme: scalar_advec = "' // &
790                        TRIM( scalar_advec ) // '"'
791       CALL message( 'check_parameters', 'PA0024', 1, 2, 0, 6, 0 )
792    ENDIF
793    IF ( scalar_advec == 'bc-scheme'  .AND.  loop_optimization == 'cache' )    &
794    THEN
795       message_string = 'advection_scheme scalar_advec = "' &
796         // TRIM( scalar_advec ) // '" not implemented for & loop_optimization = "' // &
797         TRIM( loop_optimization ) // '"'
798       CALL message( 'check_parameters', 'PA0026', 1, 2, 0, 6, 0 )
799    ENDIF
800
801    IF ( use_sgs_for_particles  .AND.  .NOT. use_upstream_for_tke .AND.        &
802         ( scalar_advec /= 'ws-scheme' .OR.                                    &
803           scalar_advec /= 'ws-scheme-mono' )                                  &
804       )  THEN
805       use_upstream_for_tke = .TRUE.
806       message_string = 'use_upstream_for_tke set .TRUE. because ' //          &
807                        'use_sgs_for_particles = .TRUE. '          //          &
808                        'and scalar_advec /= ws-scheme'
809       CALL message( 'check_parameters', 'PA0025', 0, 1, 0, 6, 0 )
810    ENDIF
811
812    IF ( use_sgs_for_particles  .AND.  curvature_solution_effects )  THEN
813       message_string = 'use_sgs_for_particles = .TRUE. not allowed with ' //  &
814                        'curvature_solution_effects = .TRUE.'
815       CALL message( 'check_parameters', 'PA0349', 1, 2, 0, 6, 0 )
816    ENDIF
817
818!
819!-- Set LOGICAL switches to enhance performance
820    IF ( momentum_advec == 'ws-scheme' )       ws_scheme_mom = .TRUE.
821    IF ( scalar_advec   == 'ws-scheme' .OR.                                   &
822         scalar_advec   == 'ws-scheme-mono' )  ws_scheme_sca = .TRUE.
823    IF ( scalar_advec   == 'ws-scheme-mono' )  monotonic_adjustment = .TRUE.
824
825
826!
827!-- Timestep schemes:
828    SELECT CASE ( TRIM( timestep_scheme ) )
829
830       CASE ( 'euler' )
831          intermediate_timestep_count_max = 1
832
833       CASE ( 'runge-kutta-2' )
834          intermediate_timestep_count_max = 2
835
836       CASE ( 'runge-kutta-3' )
837          intermediate_timestep_count_max = 3
838
839       CASE DEFAULT
840          message_string = 'unknown timestep scheme: timestep_scheme = "' //   &
841                           TRIM( timestep_scheme ) // '"'
842          CALL message( 'check_parameters', 'PA0027', 1, 2, 0, 6, 0 )
843
844    END SELECT
845
846    IF ( (momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme')   &
847         .AND. timestep_scheme(1:5) == 'runge' ) THEN
848       message_string = 'momentum advection scheme "' // &
849                        TRIM( momentum_advec ) // '" & does not work with ' // &
850                        'timestep_scheme "' // TRIM( timestep_scheme ) // '"'
851       CALL message( 'check_parameters', 'PA0029', 1, 2, 0, 6, 0 )
852    ENDIF
853
854!
855!-- Collision kernels:
856    SELECT CASE ( TRIM( collision_kernel ) )
857
858       CASE ( 'hall', 'hall_fast' )
859          hall_kernel = .TRUE.
860
861       CASE ( 'palm' )
862          palm_kernel = .TRUE.
863
864       CASE ( 'wang', 'wang_fast' )
865          wang_kernel = .TRUE.
866
867       CASE ( 'none' )
868
869
870       CASE DEFAULT
871          message_string = 'unknown collision kernel: collision_kernel = "' // &
872                           TRIM( collision_kernel ) // '"'
873          CALL message( 'check_parameters', 'PA0350', 1, 2, 0, 6, 0 )
874
875    END SELECT
876    IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
877
878    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
879         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
880!
881!--    No restart run: several initialising actions are possible
882       action = initializing_actions
883       DO WHILE ( TRIM( action ) /= '' )
884          position = INDEX( action, ' ' )
885          SELECT CASE ( action(1:position-1) )
886
887             CASE ( 'set_constant_profiles', 'set_1d-model_profiles', &
888                    'by_user', 'initialize_vortex',     'initialize_ptanom' )
889                action = action(position+1:)
890
891             CASE DEFAULT
892                message_string = 'initializing_action = "' // &
893                                 TRIM( action ) // '" unkown or not allowed'
894                CALL message( 'check_parameters', 'PA0030', 1, 2, 0, 6, 0 )
895
896          END SELECT
897       ENDDO
898    ENDIF
899
900    IF ( TRIM( initializing_actions ) == 'initialize_vortex' .AND. &
901         conserve_volume_flow ) THEN
902         message_string = 'initializing_actions = "initialize_vortex"' // &
903                        ' ist not allowed with conserve_volume_flow = .T.'
904       CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
905    ENDIF       
906
907
908    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND. &
909         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
910       message_string = 'initializing_actions = "set_constant_profiles"' // &
911                        ' and "set_1d-model_profiles" are not allowed ' //  &
912                        'simultaneously'
913       CALL message( 'check_parameters', 'PA0031', 1, 2, 0, 6, 0 )
914    ENDIF
915
916    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND. &
917         INDEX( initializing_actions, 'by_user' ) /= 0 )  THEN
918       message_string = 'initializing_actions = "set_constant_profiles"' // &
919                        ' and "by_user" are not allowed simultaneously'
920       CALL message( 'check_parameters', 'PA0032', 1, 2, 0, 6, 0 )
921    ENDIF
922
923    IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND. &
924         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
925       message_string = 'initializing_actions = "by_user" and ' // &
926                        '"set_1d-model_profiles" are not allowed simultaneously'
927       CALL message( 'check_parameters', 'PA0033', 1, 2, 0, 6, 0 )
928    ENDIF
929
930    IF ( cloud_physics  .AND.  .NOT. humidity )  THEN
931       WRITE( message_string, * ) 'cloud_physics = ', cloud_physics, ' is ', &
932              'not allowed with humidity = ', humidity
933       CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
934    ENDIF
935
936    IF ( precipitation  .AND.  .NOT.  cloud_physics )  THEN
937       WRITE( message_string, * ) 'precipitation = ', precipitation, ' is ', &
938              'not allowed with cloud_physics = ', cloud_physics
939       CALL message( 'check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
940    ENDIF
941
942    IF ( humidity  .AND.  sloping_surface )  THEN
943       message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' // &
944                        'are not allowed simultaneously'
945       CALL message( 'check_parameters', 'PA0036', 1, 2, 0, 6, 0 )
946    ENDIF
947
948    IF ( passive_scalar  .AND.  humidity )  THEN
949       message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' // &
950                        'is not allowed simultaneously'
951       CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
952    ENDIF
953
954    IF ( plant_canopy )  THEN
955   
956       IF ( canopy_drag_coeff == 0.0_wp )  THEN
957          message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// &
958                           'coefficient & given value is canopy_drag_coeff = 0.0'
959          CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
960       ENDIF
961   
962       IF ( ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad == 9999999.9_wp )  .OR.&
963              beta_lad /= 9999999.9_wp   .AND.  alpha_lad == 9999999.9_wp )  THEN
964          message_string = 'using the beta function for the construction ' //  &
965                           'of the leaf area density profile requires '    //  &
966                           'both alpha_lad and beta_lad to be /= 9999999.9'
967          CALL message( 'check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
968       ENDIF
969   
970       IF ( calc_beta_lad_profile  .AND.  lai_beta == 0.0_wp )  THEN
971          message_string = 'using the beta function for the construction ' //  &
972                           'of the leaf area density profile requires '    //  &
973                           'a non-zero lai_beta, but given value is '      //  &
974                           'lai_beta = 0.0'
975          CALL message( 'check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
976       ENDIF
977
978       IF ( calc_beta_lad_profile  .AND.  lad_surface /= 0.0_wp )  THEN
979          message_string = 'simultaneous setting of alpha_lad /= 9999999.9' // &
980                           'and lad_surface /= 0.0 is not possible, '       // &
981                           'use either vertical gradients or the beta '     // &
982                           'function for the construction of the leaf area '// &
983                           'density profile'
984          CALL message( 'check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
985       ENDIF
986
987       IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
988          message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
989                           ' seifert_beheng'
990          CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
991       ENDIF
992
993    ENDIF
994
995
996    IF ( land_surface )  THEN
997 
998!
999!--    Dirichlet boundary conditions are required as the surface fluxes are
1000!--    calculated from the temperature/humidity gradients in the land surface
1001!--    model
1002       IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' )  THEN
1003          message_string = 'lsm requires setting of'//                         &
1004                           'bc_pt_b = "dirichlet" and '//                      &
1005                           'bc_q_b  = "dirichlet"'
1006          CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
1007       ENDIF
1008
1009       IF ( .NOT. prandtl_layer )  THEN
1010          message_string = 'lsm requires '//                                   &
1011                           'prandtl_layer = .T.'
1012          CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
1013       ENDIF
1014
1015       IF ( veg_type == 0 )  THEN
1016          IF ( SUM(root_fraction) /= 1.0_wp)  THEN
1017             message_string = 'veg_type = 0 (user_defined)'//                  &
1018                              'requires setting of root_fraction(0:3)'//       &
1019                              '/= 9999999.9 and SUM(root_fraction) = 1'
1020             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1021          ENDIF
1022 
1023          IF ( min_canopy_resistance == 9999999.9_wp)  THEN
1024             message_string = 'veg_type = 0 (user defined)'//                  &
1025                              'requires setting of min_canopy_resistance'//    &
1026                              '/= 9999999.9'
1027             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1028          ENDIF
1029
1030          IF ( leaf_area_index == 9999999.9_wp)  THEN
1031             message_string = 'veg_type = 0 (user_defined)'//                  &
1032                              'requires setting of leaf_area_index'//          &
1033                              '/= 9999999.9'
1034             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1035          ENDIF
1036
1037          IF ( vegetation_coverage == 9999999.9_wp)  THEN
1038             message_string = 'veg_type = 0 (user_defined)'//                  &
1039                              'requires setting of vegetation_coverage'//      &
1040                              '/= 9999999.9'
1041             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1042          ENDIF
1043
1044          IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
1045             message_string = 'veg_type = 0 (user_defined)'//                  &
1046                              'requires setting of'//                          &
1047                              'canopy_resistance_coefficient /= 9999999.9'
1048             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1049          ENDIF
1050
1051          IF ( lambda_surface_stable == 9999999.9_wp)  THEN
1052             message_string = 'veg_type = 0 (user_defined)'//                  &
1053                              'requires setting of lambda_surface_stable'//    &
1054                              '/= 9999999.9'
1055             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1056          ENDIF
1057
1058          IF ( lambda_surface_unstable == 9999999.9_wp)  THEN
1059             message_string = 'veg_type = 0 (user_defined)'//                  &
1060                              'requires setting of lambda_surface_unstable'//  &
1061                              '/= 9999999.9'
1062             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1063          ENDIF
1064
1065          IF ( f_shortwave_incoming == 9999999.9_wp)  THEN
1066             message_string = 'veg_type = 0 (user_defined)'//                  &
1067                              'requires setting of f_shortwave_incoming'//     &
1068                              '/= 9999999.9'
1069             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1070          ENDIF
1071
1072          IF ( z0_eb == 9999999.9_wp)  THEN
1073             message_string = 'veg_type = 0 (user_defined)'//                  &
1074                              'requires setting of z0_eb'//                   &
1075                              '/= 9999999.9'
1076             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1077          ENDIF
1078
1079          IF ( z0h_eb == 9999999.9_wp)  THEN
1080             message_string = 'veg_type = 0 (user_defined)'//                  &
1081                              'requires setting of z0h_eb'//                  &
1082                              '/= 9999999.9'
1083             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1084          ENDIF
1085
1086
1087       ENDIF
1088
1089       IF ( soil_type == 0 )  THEN
1090
1091          IF ( alpha_vangenuchten == 9999999.9_wp)  THEN
1092             message_string = 'soil_type = 0 (user_defined)'//                 &
1093                              'requires setting of alpha_vangenuchten'//       &
1094                              '/= 9999999.9'
1095             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1096          ENDIF
1097
1098          IF ( l_vangenuchten == 9999999.9_wp)  THEN
1099             message_string = 'soil_type = 0 (user_defined)'//                 &
1100                              'requires setting of l_vangenuchten'//           &
1101                              '/= 9999999.9'
1102             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1103          ENDIF
1104
1105          IF ( n_vangenuchten == 9999999.9_wp)  THEN
1106             message_string = 'soil_type = 0 (user_defined)'//                 &
1107                              'requires setting of n_vangenuchten'//           &
1108                              '/= 9999999.9'
1109             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1110          ENDIF
1111
1112          IF ( hydraulic_conductivity == 9999999.9_wp)  THEN
1113             message_string = 'soil_type = 0 (user_defined)'//                 &
1114                              'requires setting of hydraulic_conductivity'//   &
1115                              '/= 9999999.9'
1116             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1117          ENDIF
1118
1119          IF ( saturation_moisture == 9999999.9_wp)  THEN
1120             message_string = 'soil_type = 0 (user_defined)'//                 &
1121                              'requires setting of saturation_moisture'//      &
1122                              '/= 9999999.9'
1123             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1124          ENDIF
1125
1126          IF ( field_capacity == 9999999.9_wp)  THEN
1127             message_string = 'soil_type = 0 (user_defined)'//                 &
1128                              'requires setting of field_capacity'//           &
1129                              '/= 9999999.9'
1130             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1131          ENDIF
1132
1133          IF ( wilting_point == 9999999.9_wp)  THEN
1134             message_string = 'soil_type = 0 (user_defined)'//                 &
1135                              'requires setting of wilting_point'//            &
1136                              '/= 9999999.9'
1137             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1138          ENDIF
1139
1140          IF ( residual_moisture == 9999999.9_wp)  THEN
1141             message_string = 'soil_type = 0 (user_defined)'//                 &
1142                              'requires setting of residual_moisture'//        &
1143                              '/= 9999999.9'
1144             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1145          ENDIF
1146
1147       ENDIF
1148
1149       IF ( .NOT. radiation )  THEN
1150          message_string = 'lsm requires '//                                   &
1151                           'radiation = .T.'
1152          CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
1153       ENDIF
1154
1155    END IF
1156
1157    IF ( radiation )  THEN
1158       IF ( radiation_scheme /= 'constant'  .AND.                              &
1159            radiation_scheme /= 'clear-sky' .AND.                              &
1160            radiation_scheme /= 'rrtmg' )  THEN
1161          message_string = 'unknown radiation_scheme = '//                     &
1162                           TRIM( radiation_scheme )
1163          CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
1164       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
1165#if ! defined ( __rrtmg )
1166          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
1167                           'compilation of PALM with pre-processor ' //        &
1168                           'directive -D__rrtmg'
1169          CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
1170#endif
1171#if defined ( __rrtmg ) && ! defined( __netcdf )
1172          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
1173                           'the use of NetCDF (preprocessor directive ' //     &
1174                           '-D__netcdf'
1175          CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
1176#endif
1177
1178       ENDIF
1179       IF ( albedo_type == 0 .AND. albedo == 9999999.9_wp .AND.                &
1180            radiation_scheme == 'clear-sky')  THEN
1181          message_string = 'radiation_scheme = "clear-sky" in combination' //  & 
1182                           'with albedo_type = 0 requires setting of albedo'// &
1183                           ' /= 9999999.9'
1184          CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
1185       ENDIF
1186       IF ( albedo_type == 0 .AND. radiation_scheme == 'rrtmg' .AND.           &
1187          (    albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
1188          .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& 
1189          ) ) THEN
1190          message_string = 'radiation_scheme = "rrtmg" in combination' //      & 
1191                           'with albedo_type = 0 requires setting of ' //      &
1192                           'albedo_lw_dif /= 9999999.9' //                     &
1193                           'albedo_lw_dir /= 9999999.9' //                     &
1194                           'albedo_sw_dif /= 9999999.9 and' //                 &
1195                           'albedo_sw_dir /= 9999999.9'
1196          CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
1197       ENDIF
1198
1199    ENDIF
1200
1201
1202    IF ( .NOT. ( loop_optimization == 'cache'  .OR.                            &
1203                 loop_optimization == 'vector' )                               &
1204         .AND.  cloud_physics  .AND.  icloud_scheme == 0 )  THEN
1205       message_string = 'cloud_scheme = seifert_beheng requires ' // &
1206                        'loop_optimization = "cache" or "vector"'
1207       CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
1208    ENDIF 
1209
1210!
1211!-- In case of no model continuation run, check initialising parameters and
1212!-- deduce further quantities
1213    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1214
1215!
1216!--    Initial profiles for 1D and 3D model, respectively (u,v further below)
1217       pt_init = pt_surface
1218       IF ( humidity )  THEN
1219          q_init  = q_surface
1220       ENDIF
1221       IF ( ocean )           sa_init = sa_surface
1222       IF ( passive_scalar )  q_init  = s_surface
1223
1224!
1225!--
1226!--    If required, compute initial profile of the geostrophic wind
1227!--    (component ug)
1228       i = 1
1229       gradient = 0.0_wp
1230
1231       IF ( .NOT. ocean )  THEN
1232
1233          ug_vertical_gradient_level_ind(1) = 0
1234          ug(0) = ug_surface
1235          DO  k = 1, nzt+1
1236             IF ( i < 11 ) THEN
1237                IF ( ug_vertical_gradient_level(i) < zu(k)  .AND. &
1238                     ug_vertical_gradient_level(i) >= 0.0_wp )  THEN
1239                   gradient = ug_vertical_gradient(i) / 100.0_wp
1240                   ug_vertical_gradient_level_ind(i) = k - 1
1241                   i = i + 1
1242                ENDIF
1243             ENDIF       
1244             IF ( gradient /= 0.0_wp )  THEN
1245                IF ( k /= 1 )  THEN
1246                   ug(k) = ug(k-1) + dzu(k) * gradient
1247                ELSE
1248                   ug(k) = ug_surface + dzu(k) * gradient
1249                ENDIF
1250             ELSE
1251                ug(k) = ug(k-1)
1252             ENDIF
1253          ENDDO
1254
1255       ELSE
1256
1257          ug_vertical_gradient_level_ind(1) = nzt+1
1258          ug(nzt+1) = ug_surface
1259          DO  k = nzt, nzb, -1
1260             IF ( i < 11 ) THEN
1261                IF ( ug_vertical_gradient_level(i) > zu(k)  .AND. &
1262                     ug_vertical_gradient_level(i) <= 0.0_wp )  THEN
1263                   gradient = ug_vertical_gradient(i) / 100.0_wp
1264                   ug_vertical_gradient_level_ind(i) = k + 1
1265                   i = i + 1
1266                ENDIF
1267             ENDIF
1268             IF ( gradient /= 0.0_wp )  THEN
1269                IF ( k /= nzt )  THEN
1270                   ug(k) = ug(k+1) - dzu(k+1) * gradient
1271                ELSE
1272                   ug(k)   = ug_surface - 0.5_wp * dzu(k+1) * gradient
1273                   ug(k+1) = ug_surface + 0.5_wp * dzu(k+1) * gradient
1274                ENDIF
1275             ELSE
1276                ug(k) = ug(k+1)
1277             ENDIF
1278          ENDDO
1279
1280       ENDIF
1281
1282!
1283!--    In case of no given gradients for ug, choose a zero gradient
1284       IF ( ug_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1285          ug_vertical_gradient_level(1) = 0.0_wp
1286       ENDIF 
1287
1288!
1289!--
1290!--    If required, compute initial profile of the geostrophic wind
1291!--    (component vg)
1292       i = 1
1293       gradient = 0.0_wp
1294
1295       IF ( .NOT. ocean )  THEN
1296
1297          vg_vertical_gradient_level_ind(1) = 0
1298          vg(0) = vg_surface
1299          DO  k = 1, nzt+1
1300             IF ( i < 11 ) THEN
1301                IF ( vg_vertical_gradient_level(i) < zu(k)  .AND. &
1302                     vg_vertical_gradient_level(i) >= 0.0_wp )  THEN
1303                   gradient = vg_vertical_gradient(i) / 100.0_wp
1304                   vg_vertical_gradient_level_ind(i) = k - 1
1305                   i = i + 1
1306                ENDIF
1307             ENDIF
1308             IF ( gradient /= 0.0_wp )  THEN
1309                IF ( k /= 1 )  THEN
1310                   vg(k) = vg(k-1) + dzu(k) * gradient
1311                ELSE
1312                   vg(k) = vg_surface + dzu(k) * gradient
1313                ENDIF
1314             ELSE
1315                vg(k) = vg(k-1)
1316             ENDIF
1317          ENDDO
1318
1319       ELSE
1320
1321          vg_vertical_gradient_level_ind(1) = nzt+1
1322          vg(nzt+1) = vg_surface
1323          DO  k = nzt, nzb, -1
1324             IF ( i < 11 ) THEN
1325                IF ( vg_vertical_gradient_level(i) > zu(k)  .AND. &
1326                     vg_vertical_gradient_level(i) <= 0.0_wp )  THEN
1327                   gradient = vg_vertical_gradient(i) / 100.0_wp
1328                   vg_vertical_gradient_level_ind(i) = k + 1
1329                   i = i + 1
1330                ENDIF
1331             ENDIF
1332             IF ( gradient /= 0.0_wp )  THEN
1333                IF ( k /= nzt )  THEN
1334                   vg(k) = vg(k+1) - dzu(k+1) * gradient
1335                ELSE
1336                   vg(k)   = vg_surface - 0.5_wp * dzu(k+1) * gradient
1337                   vg(k+1) = vg_surface + 0.5_wp * dzu(k+1) * gradient
1338                ENDIF
1339             ELSE
1340                vg(k) = vg(k+1)
1341             ENDIF
1342          ENDDO
1343
1344       ENDIF
1345
1346!
1347!--    In case of no given gradients for vg, choose a zero gradient
1348       IF ( vg_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1349          vg_vertical_gradient_level(1) = 0.0_wp
1350       ENDIF
1351
1352!
1353!--    Let the initial wind profiles be the calculated ug/vg profiles or
1354!--    interpolate them from wind profile data (if given)
1355       IF ( u_profile(1) == 9999999.9_wp  .AND.  v_profile(1) == 9999999.9_wp )  THEN
1356
1357          u_init = ug
1358          v_init = vg
1359
1360       ELSEIF ( u_profile(1) == 0.0_wp  .AND.  v_profile(1) == 0.0_wp )  THEN
1361
1362          IF ( uv_heights(1) /= 0.0_wp )  THEN
1363             message_string = 'uv_heights(1) must be 0.0'
1364             CALL message( 'check_parameters', 'PA0345', 1, 2, 0, 6, 0 )
1365          ENDIF
1366
1367          use_prescribed_profile_data = .TRUE.
1368
1369          kk = 1
1370          u_init(0) = 0.0_wp
1371          v_init(0) = 0.0_wp
1372
1373          DO  k = 1, nz+1
1374
1375             IF ( kk < 100 )  THEN
1376                DO WHILE ( uv_heights(kk+1) <= zu(k) )
1377                   kk = kk + 1
1378                   IF ( kk == 100 )  EXIT
1379                ENDDO
1380             ENDIF
1381
1382             IF ( kk < 100 .AND. uv_heights(kk+1) /= 9999999.9_wp )  THEN
1383                u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1384                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1385                                       ( u_profile(kk+1) - u_profile(kk) )
1386                v_init(k) = v_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1387                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1388                                       ( v_profile(kk+1) - v_profile(kk) )
1389             ELSE
1390                u_init(k) = u_profile(kk)
1391                v_init(k) = v_profile(kk)
1392             ENDIF
1393
1394          ENDDO
1395
1396       ELSE
1397
1398          message_string = 'u_profile(1) and v_profile(1) must be 0.0'
1399          CALL message( 'check_parameters', 'PA0346', 1, 2, 0, 6, 0 )
1400
1401       ENDIF
1402
1403!
1404!--    Compute initial temperature profile using the given temperature gradients
1405       IF ( .NOT. neutral )  THEN
1406
1407          i = 1
1408          gradient = 0.0_wp
1409
1410          IF ( .NOT. ocean )  THEN
1411
1412             pt_vertical_gradient_level_ind(1) = 0
1413             DO  k = 1, nzt+1
1414                IF ( i < 11 ) THEN
1415                   IF ( pt_vertical_gradient_level(i) < zu(k)  .AND. &
1416                        pt_vertical_gradient_level(i) >= 0.0_wp )  THEN
1417                      gradient = pt_vertical_gradient(i) / 100.0_wp
1418                      pt_vertical_gradient_level_ind(i) = k - 1
1419                      i = i + 1
1420                   ENDIF
1421                ENDIF
1422                IF ( gradient /= 0.0_wp )  THEN
1423                   IF ( k /= 1 )  THEN
1424                      pt_init(k) = pt_init(k-1) + dzu(k) * gradient
1425                   ELSE
1426                      pt_init(k) = pt_surface   + dzu(k) * gradient
1427                   ENDIF
1428                ELSE
1429                   pt_init(k) = pt_init(k-1)
1430                ENDIF
1431             ENDDO
1432
1433          ELSE
1434
1435             pt_vertical_gradient_level_ind(1) = nzt+1
1436             DO  k = nzt, 0, -1
1437                IF ( i < 11 ) THEN
1438                   IF ( pt_vertical_gradient_level(i) > zu(k)  .AND. &
1439                        pt_vertical_gradient_level(i) <= 0.0_wp )  THEN
1440                      gradient = pt_vertical_gradient(i) / 100.0_wp
1441                      pt_vertical_gradient_level_ind(i) = k + 1
1442                      i = i + 1
1443                   ENDIF
1444                ENDIF
1445                IF ( gradient /= 0.0_wp )  THEN
1446                   IF ( k /= nzt )  THEN
1447                      pt_init(k) = pt_init(k+1) - dzu(k+1) * gradient
1448                   ELSE
1449                      pt_init(k)   = pt_surface - 0.5_wp * dzu(k+1) * gradient
1450                      pt_init(k+1) = pt_surface + 0.5_wp * dzu(k+1) * gradient
1451                   ENDIF
1452                ELSE
1453                   pt_init(k) = pt_init(k+1)
1454                ENDIF
1455             ENDDO
1456
1457          ENDIF
1458
1459       ENDIF
1460
1461!
1462!--    In case of no given temperature gradients, choose gradient of neutral
1463!--    stratification
1464       IF ( pt_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1465          pt_vertical_gradient_level(1) = 0.0_wp
1466       ENDIF
1467
1468!
1469!--    Store temperature gradient at the top boundary for possible Neumann
1470!--    boundary condition
1471       bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
1472
1473!
1474!--    If required, compute initial humidity or scalar profile using the given
1475!--    humidity/scalar gradient. In case of scalar transport, initially store
1476!--    values of the scalar parameters on humidity parameters
1477       IF ( passive_scalar )  THEN
1478          bc_q_b                    = bc_s_b
1479          bc_q_t                    = bc_s_t
1480          q_surface                 = s_surface
1481          q_surface_initial_change  = s_surface_initial_change
1482          q_vertical_gradient       = s_vertical_gradient
1483          q_vertical_gradient_level = s_vertical_gradient_level
1484          surface_waterflux         = surface_scalarflux
1485          wall_humidityflux         = wall_scalarflux
1486       ENDIF
1487
1488       IF ( humidity  .OR.  passive_scalar )  THEN
1489
1490          i = 1
1491          gradient = 0.0_wp
1492          q_vertical_gradient_level_ind(1) = 0
1493          DO  k = 1, nzt+1
1494             IF ( i < 11 ) THEN
1495                IF ( q_vertical_gradient_level(i) < zu(k)  .AND. &
1496                     q_vertical_gradient_level(i) >= 0.0_wp )  THEN
1497                   gradient = q_vertical_gradient(i) / 100.0_wp
1498                   q_vertical_gradient_level_ind(i) = k - 1
1499                   i = i + 1
1500                ENDIF
1501             ENDIF
1502             IF ( gradient /= 0.0_wp )  THEN
1503                IF ( k /= 1 )  THEN
1504                   q_init(k) = q_init(k-1) + dzu(k) * gradient
1505                ELSE
1506                   q_init(k) = q_init(k-1) + dzu(k) * gradient
1507                ENDIF
1508             ELSE
1509                q_init(k) = q_init(k-1)
1510             ENDIF
1511!
1512!--          Avoid negative humidities
1513             IF ( q_init(k) < 0.0_wp )  THEN
1514                q_init(k) = 0.0_wp
1515             ENDIF
1516          ENDDO
1517
1518!
1519!--       In case of no given humidity gradients, choose zero gradient
1520!--       conditions
1521          IF ( q_vertical_gradient_level(1) == -1.0_wp )  THEN
1522             q_vertical_gradient_level(1) = 0.0_wp
1523          ENDIF
1524!
1525!--       Store humidity, rain water content and rain drop concentration
1526!--       gradient at the top boundary for possile Neumann boundary condition
1527          bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
1528       ENDIF
1529
1530!
1531!--    If required, compute initial salinity profile using the given salinity
1532!--    gradients
1533       IF ( ocean )  THEN
1534
1535          i = 1
1536          gradient = 0.0_wp
1537
1538          sa_vertical_gradient_level_ind(1) = nzt+1
1539          DO  k = nzt, 0, -1
1540             IF ( i < 11 ) THEN
1541                IF ( sa_vertical_gradient_level(i) > zu(k)  .AND. &
1542                     sa_vertical_gradient_level(i) <= 0.0_wp )  THEN
1543                   gradient = sa_vertical_gradient(i) / 100.0_wp
1544                   sa_vertical_gradient_level_ind(i) = k + 1
1545                   i = i + 1
1546                ENDIF
1547             ENDIF
1548             IF ( gradient /= 0.0_wp )  THEN
1549                IF ( k /= nzt )  THEN
1550                   sa_init(k) = sa_init(k+1) - dzu(k+1) * gradient
1551                ELSE
1552                   sa_init(k)   = sa_surface - 0.5_wp * dzu(k+1) * gradient
1553                   sa_init(k+1) = sa_surface + 0.5_wp * dzu(k+1) * gradient
1554                ENDIF
1555             ELSE
1556                sa_init(k) = sa_init(k+1)
1557             ENDIF
1558          ENDDO
1559
1560       ENDIF
1561
1562         
1563    ENDIF
1564
1565!
1566!-- Check if the control parameter use_subsidence_tendencies is used correctly
1567    IF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_subsidence )  THEN
1568       message_string = 'The usage of use_subsidence_tendencies ' // &
1569                            'requires large_scale_subsidence = .T..'
1570       CALL message( 'check_parameters', 'PA0396', 1, 2, 0, 6, 0 )
1571    ELSEIF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_forcing )  THEN
1572       message_string = 'The usage of use_subsidence_tendencies ' // &
1573                            'requires large_scale_forcing = .T..'
1574       CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 )
1575    ENDIF
1576
1577!
1578!-- Initialize large scale subsidence if required
1579    If ( large_scale_subsidence )  THEN
1580       IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp .AND. &
1581                                     .NOT. large_scale_forcing )  THEN
1582          CALL init_w_subsidence
1583       ENDIF
1584!
1585!--    In case large_scale_forcing is used, profiles for subsidence velocity
1586!--    are read in from file LSF_DATA
1587
1588       IF ( subs_vertical_gradient_level(1) == -9999999.9_wp .AND. &
1589                                     .NOT. large_scale_forcing )  THEN
1590          message_string = 'There is no default large scale vertical ' // &
1591                           'velocity profile set. Specify the subsidence ' // &
1592                           'velocity profile via subs_vertical_gradient and ' // &
1593                           'subs_vertical_gradient_level.'
1594          CALL message( 'check_parameters', 'PA0380', 1, 2, 0, 6, 0 )
1595       ENDIF
1596    ELSE
1597        IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp )  THEN
1598           message_string = 'Enable usage of large scale subsidence by ' // &
1599                            'setting large_scale_subsidence = .T..'
1600          CALL message( 'check_parameters', 'PA0381', 1, 2, 0, 6, 0 )
1601        ENDIF
1602    ENDIF   
1603
1604!
1605!-- Compute Coriolis parameter
1606    f  = 2.0_wp * omega * SIN( phi / 180.0_wp * pi )
1607    fs = 2.0_wp * omega * COS( phi / 180.0_wp * pi )
1608
1609!
1610!-- Check and set buoyancy related parameters and switches
1611    IF ( reference_state == 'horizontal_average' )  THEN
1612       CONTINUE
1613    ELSEIF ( reference_state == 'initial_profile' )  THEN
1614       use_initial_profile_as_reference = .TRUE.
1615    ELSEIF ( reference_state == 'single_value' )  THEN
1616       use_single_reference_value = .TRUE.
1617       IF ( pt_reference == 9999999.9_wp )  pt_reference = pt_surface
1618       vpt_reference = pt_reference * ( 1.0_wp + 0.61_wp * q_surface )
1619    ELSE
1620       message_string = 'illegal value for reference_state: "' // &
1621                        TRIM( reference_state ) // '"'
1622       CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
1623    ENDIF
1624
1625!
1626!-- Ocean runs always use reference values in the buoyancy term
1627    IF ( ocean )  THEN
1628       reference_state = 'single_value'
1629       use_single_reference_value = .TRUE.
1630    ENDIF
1631
1632!
1633!-- Sign of buoyancy/stability terms
1634    IF ( ocean )  atmos_ocean_sign = -1.0_wp
1635
1636!
1637!-- Ocean version must use flux boundary conditions at the top
1638    IF ( ocean .AND. .NOT. use_top_fluxes )  THEN
1639       message_string = 'use_top_fluxes must be .TRUE. in ocean mode'
1640       CALL message( 'check_parameters', 'PA0042', 1, 2, 0, 6, 0 )
1641    ENDIF
1642
1643!
1644!-- In case of a given slope, compute the relevant quantities
1645    IF ( alpha_surface /= 0.0_wp )  THEN
1646       IF ( ABS( alpha_surface ) > 90.0_wp )  THEN
1647          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface, &
1648                                     ' ) must be < 90.0'
1649          CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 )
1650       ENDIF
1651       sloping_surface = .TRUE.
1652       cos_alpha_surface = COS( alpha_surface / 180.0_wp * pi )
1653       sin_alpha_surface = SIN( alpha_surface / 180.0_wp * pi )
1654    ENDIF
1655
1656!
1657!-- Check time step and cfl_factor
1658    IF ( dt /= -1.0_wp )  THEN
1659       IF ( dt <= 0.0_wp  .AND.  dt /= -1.0_wp )  THEN
1660          WRITE( message_string, * ) 'dt = ', dt , ' <= 0.0'
1661          CALL message( 'check_parameters', 'PA0044', 1, 2, 0, 6, 0 )
1662       ENDIF
1663       dt_3d = dt
1664       dt_fixed = .TRUE.
1665    ENDIF
1666
1667    IF ( cfl_factor <= 0.0_wp  .OR.  cfl_factor > 1.0_wp )  THEN
1668       IF ( cfl_factor == -1.0_wp )  THEN
1669          IF ( timestep_scheme == 'runge-kutta-2' )  THEN
1670             cfl_factor = 0.8_wp
1671          ELSEIF ( timestep_scheme == 'runge-kutta-3' )  THEN
1672             cfl_factor = 0.9_wp
1673          ELSE
1674             cfl_factor = 0.9_wp
1675          ENDIF
1676       ELSE
1677          WRITE( message_string, * ) 'cfl_factor = ', cfl_factor, &
1678                 ' out of range & 0.0 < cfl_factor <= 1.0 is required'
1679          CALL message( 'check_parameters', 'PA0045', 1, 2, 0, 6, 0 )
1680       ENDIF
1681    ENDIF
1682
1683!
1684!-- Store simulated time at begin
1685    simulated_time_at_begin = simulated_time
1686
1687!
1688!-- Store reference time for coupled runs and change the coupling flag,
1689!-- if ...
1690    IF ( simulated_time == 0.0_wp )  THEN
1691       IF ( coupling_start_time == 0.0_wp )  THEN
1692          time_since_reference_point = 0.0_wp
1693       ELSEIF ( time_since_reference_point < 0.0_wp )  THEN
1694          run_coupled = .FALSE.
1695       ENDIF
1696    ENDIF
1697
1698!
1699!-- Set wind speed in the Galilei-transformed system
1700    IF ( galilei_transformation )  THEN
1701       IF ( use_ug_for_galilei_tr .AND.                     &
1702            ug_vertical_gradient_level(1) == 0.0_wp  .AND.  &
1703            ug_vertical_gradient(1) == 0.0_wp  .AND.        & 
1704            vg_vertical_gradient_level(1) == 0.0_wp  .AND.  &
1705            vg_vertical_gradient(1) == 0.0_wp )  THEN
1706          u_gtrans = ug_surface * 0.6_wp
1707          v_gtrans = vg_surface * 0.6_wp
1708       ELSEIF ( use_ug_for_galilei_tr  .AND.                     &
1709                ( ug_vertical_gradient_level(1) /= 0.0_wp  .OR.  &
1710                ug_vertical_gradient(1) /= 0.0_wp ) )  THEN
1711          message_string = 'baroclinity (ug) not allowed simultaneously' // &
1712                           ' with galilei transformation'
1713          CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
1714       ELSEIF ( use_ug_for_galilei_tr  .AND.                     &
1715                ( vg_vertical_gradient_level(1) /= 0.0_wp  .OR.  &
1716                vg_vertical_gradient(1) /= 0.0_wp ) )  THEN
1717          message_string = 'baroclinity (vg) not allowed simultaneously' // &
1718                           ' with galilei transformation'
1719          CALL message( 'check_parameters', 'PA0047', 1, 2, 0, 6, 0 )
1720       ELSE
1721          message_string = 'variable translation speed used for galilei-' // &
1722             'transformation, which may cause & instabilities in stably ' // &
1723             'stratified regions'
1724          CALL message( 'check_parameters', 'PA0048', 0, 1, 0, 6, 0 )
1725       ENDIF
1726    ENDIF
1727
1728!
1729!-- In case of using a prandtl-layer, calculated (or prescribed) surface
1730!-- fluxes have to be used in the diffusion-terms
1731    IF ( prandtl_layer )  use_surface_fluxes = .TRUE.
1732
1733!
1734!-- Check boundary conditions and set internal variables:
1735!-- Lateral boundary conditions
1736    IF ( bc_lr /= 'cyclic'  .AND.  bc_lr /= 'dirichlet/radiation'  .AND. &
1737         bc_lr /= 'radiation/dirichlet' )  THEN
1738       message_string = 'unknown boundary condition: bc_lr = "' // &
1739                        TRIM( bc_lr ) // '"'
1740       CALL message( 'check_parameters', 'PA0049', 1, 2, 0, 6, 0 )
1741    ENDIF
1742    IF ( bc_ns /= 'cyclic'  .AND.  bc_ns /= 'dirichlet/radiation'  .AND. &
1743         bc_ns /= 'radiation/dirichlet' )  THEN
1744       message_string = 'unknown boundary condition: bc_ns = "' // &
1745                        TRIM( bc_ns ) // '"'
1746       CALL message( 'check_parameters', 'PA0050', 1, 2, 0, 6, 0 )
1747    ENDIF
1748
1749!
1750!-- Internal variables used for speed optimization in if clauses
1751    IF ( bc_lr /= 'cyclic' )               bc_lr_cyc    = .FALSE.
1752    IF ( bc_lr == 'dirichlet/radiation' )  bc_lr_dirrad = .TRUE.
1753    IF ( bc_lr == 'radiation/dirichlet' )  bc_lr_raddir = .TRUE.
1754    IF ( bc_ns /= 'cyclic' )               bc_ns_cyc    = .FALSE.
1755    IF ( bc_ns == 'dirichlet/radiation' )  bc_ns_dirrad = .TRUE.
1756    IF ( bc_ns == 'radiation/dirichlet' )  bc_ns_raddir = .TRUE.
1757
1758!
1759!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
1760!-- Willimas or Wicker - Skamarock advection scheme. Several schemes
1761!-- and tools do not work with non-cyclic boundary conditions.
1762    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
1763       IF ( psolver(1:9) /= 'multigrid' )  THEN
1764          message_string = 'non-cyclic lateral boundaries do not allow ' // &
1765                           'psolver = "' // TRIM( psolver ) // '"'
1766          CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
1767       ENDIF
1768       IF ( momentum_advec /= 'pw-scheme' .AND.                               &
1769          ( momentum_advec /= 'ws-scheme' .AND.                               &
1770            momentum_advec /= 'ws-scheme-mono' )                              &
1771          )  THEN
1772
1773          message_string = 'non-cyclic lateral boundaries do not allow ' // &
1774                           'momentum_advec = "' // TRIM( momentum_advec ) // '"'
1775          CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
1776       ENDIF
1777       IF ( scalar_advec /= 'pw-scheme' .AND.                                  &
1778          ( scalar_advec /= 'ws-scheme' .AND.                                  &
1779            scalar_advec /= 'ws-scheme-mono' )                                 &
1780          )  THEN
1781          message_string = 'non-cyclic lateral boundaries do not allow ' // &
1782                           'scalar_advec = "' // TRIM( scalar_advec ) // '"'
1783          CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
1784       ENDIF
1785       IF ( galilei_transformation )  THEN
1786          message_string = 'non-cyclic lateral boundaries do not allow ' // &
1787                           'galilei_transformation = .T.'
1788          CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 )
1789       ENDIF
1790    ENDIF
1791
1792!
1793!-- Bottom boundary condition for the turbulent Kinetic energy
1794    IF ( bc_e_b == 'neumann' )  THEN
1795       ibc_e_b = 1
1796    ELSEIF ( bc_e_b == '(u*)**2+neumann' )  THEN
1797       ibc_e_b = 2
1798       IF ( .NOT. prandtl_layer )  THEN
1799          bc_e_b = 'neumann'
1800          ibc_e_b = 1
1801          message_string = 'boundary condition bc_e_b changed to "' // &
1802                           TRIM( bc_e_b ) // '"'
1803          CALL message( 'check_parameters', 'PA0057', 0, 1, 0, 6, 0 )
1804       ENDIF
1805    ELSE
1806       message_string = 'unknown boundary condition: bc_e_b = "' // &
1807                        TRIM( bc_e_b ) // '"'
1808       CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 )
1809    ENDIF
1810
1811!
1812!-- Boundary conditions for perturbation pressure
1813    IF ( bc_p_b == 'dirichlet' )  THEN
1814       ibc_p_b = 0
1815    ELSEIF ( bc_p_b == 'neumann' )  THEN
1816       ibc_p_b = 1
1817    ELSE
1818       message_string = 'unknown boundary condition: bc_p_b = "' // &
1819                        TRIM( bc_p_b ) // '"'
1820       CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
1821    ENDIF
1822
1823    IF ( bc_p_t == 'dirichlet' )  THEN
1824       ibc_p_t = 0
1825    ELSEIF ( bc_p_t == 'neumann' )  THEN
1826       ibc_p_t = 1
1827    ELSE
1828       message_string = 'unknown boundary condition: bc_p_t = "' // &
1829                        TRIM( bc_p_t ) // '"'
1830       CALL message( 'check_parameters', 'PA0061', 1, 2, 0, 6, 0 )
1831    ENDIF
1832
1833!
1834!-- Boundary conditions for potential temperature
1835    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1836       ibc_pt_b = 2
1837    ELSE
1838       IF ( bc_pt_b == 'dirichlet' )  THEN
1839          ibc_pt_b = 0
1840       ELSEIF ( bc_pt_b == 'neumann' )  THEN
1841          ibc_pt_b = 1
1842       ELSE
1843          message_string = 'unknown boundary condition: bc_pt_b = "' // &
1844                           TRIM( bc_pt_b ) // '"'
1845          CALL message( 'check_parameters', 'PA0062', 1, 2, 0, 6, 0 )
1846       ENDIF
1847    ENDIF
1848
1849    IF ( bc_pt_t == 'dirichlet' )  THEN
1850       ibc_pt_t = 0
1851    ELSEIF ( bc_pt_t == 'neumann' )  THEN
1852       ibc_pt_t = 1
1853    ELSEIF ( bc_pt_t == 'initial_gradient' )  THEN
1854       ibc_pt_t = 2
1855    ELSE
1856       message_string = 'unknown boundary condition: bc_pt_t = "' // &
1857                        TRIM( bc_pt_t ) // '"'
1858       CALL message( 'check_parameters', 'PA0063', 1, 2, 0, 6, 0 )
1859    ENDIF
1860
1861    IF ( surface_heatflux == 9999999.9_wp  )  THEN
1862       constant_heatflux = .FALSE.
1863       IF ( large_scale_forcing  .OR.  land_surface )  THEN
1864          IF ( ibc_pt_b == 0 )  THEN
1865             constant_heatflux = .FALSE.
1866          ELSEIF ( ibc_pt_b == 1 )  THEN
1867             constant_heatflux = .TRUE.
1868             IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND.    &
1869                  .NOT. land_surface )  THEN
1870                surface_heatflux = shf_surf(1)
1871             ELSE
1872                surface_heatflux = 0.0_wp
1873             ENDIF
1874          ENDIF
1875       ENDIF
1876    ELSE
1877        constant_heatflux = .TRUE.
1878        IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND.         &
1879             large_scale_forcing .AND. .NOT. land_surface )  THEN
1880           surface_heatflux = shf_surf(1)
1881        ENDIF
1882    ENDIF
1883
1884    IF ( top_heatflux     == 9999999.9_wp )  constant_top_heatflux = .FALSE.
1885
1886    IF ( neutral )  THEN
1887
1888       IF ( surface_heatflux /= 0.0_wp  .AND. surface_heatflux /= 9999999.9_wp ) &
1889       THEN
1890          message_string = 'heatflux must not be set for pure neutral flow'
1891          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1892       ENDIF
1893
1894       IF ( top_heatflux /= 0.0_wp  .AND.  top_heatflux /= 9999999.9_wp ) &
1895       THEN
1896          message_string = 'heatflux must not be set for pure neutral flow'
1897          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1898       ENDIF
1899
1900    ENDIF
1901
1902    IF ( top_momentumflux_u /= 9999999.9_wp  .AND.  &
1903         top_momentumflux_v /= 9999999.9_wp )  THEN
1904       constant_top_momentumflux = .TRUE.
1905    ELSEIF (  .NOT. ( top_momentumflux_u == 9999999.9_wp  .AND.  &
1906           top_momentumflux_v == 9999999.9_wp ) )  THEN
1907       message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' // &
1908                        'must be set'
1909       CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 )
1910    ENDIF
1911
1912!
1913!-- A given surface temperature implies Dirichlet boundary condition for
1914!-- temperature. In this case specification of a constant heat flux is
1915!-- forbidden.
1916    IF ( ibc_pt_b == 0  .AND.   constant_heatflux  .AND. &
1917         surface_heatflux /= 0.0_wp )  THEN
1918       message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //&
1919                        '& is not allowed with constant_heatflux = .TRUE.'
1920       CALL message( 'check_parameters', 'PA0065', 1, 2, 0, 6, 0 )
1921    ENDIF
1922    IF ( constant_heatflux  .AND.  pt_surface_initial_change /= 0.0_wp )  THEN
1923       WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo', &
1924               'wed with pt_surface_initial_change (/=0) = ', &
1925               pt_surface_initial_change
1926       CALL message( 'check_parameters', 'PA0066', 1, 2, 0, 6, 0 )
1927    ENDIF
1928
1929!
1930!-- A given temperature at the top implies Dirichlet boundary condition for
1931!-- temperature. In this case specification of a constant heat flux is
1932!-- forbidden.
1933    IF ( ibc_pt_t == 0  .AND.   constant_top_heatflux  .AND. &
1934         top_heatflux /= 0.0_wp )  THEN
1935       message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //&
1936                        '" is not allowed with constant_top_heatflux = .TRUE.'
1937       CALL message( 'check_parameters', 'PA0067', 1, 2, 0, 6, 0 )
1938    ENDIF
1939
1940!
1941!-- Boundary conditions for salinity
1942    IF ( ocean )  THEN
1943       IF ( bc_sa_t == 'dirichlet' )  THEN
1944          ibc_sa_t = 0
1945       ELSEIF ( bc_sa_t == 'neumann' )  THEN
1946          ibc_sa_t = 1
1947       ELSE
1948          message_string = 'unknown boundary condition: bc_sa_t = "' // &
1949                           TRIM( bc_sa_t ) // '"'
1950          CALL message( 'check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
1951       ENDIF
1952
1953       IF ( top_salinityflux == 9999999.9_wp )  constant_top_salinityflux = .FALSE.
1954       IF ( ibc_sa_t == 1  .AND.   top_salinityflux == 9999999.9_wp )  THEN
1955          message_string = 'boundary condition: bc_sa_t = "' // &
1956                           TRIM( bc_sa_t ) // '" requires to set ' // &
1957                           'top_salinityflux'
1958          CALL message( 'check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
1959       ENDIF
1960
1961!
1962!--    A fixed salinity at the top implies Dirichlet boundary condition for
1963!--    salinity. In this case specification of a constant salinity flux is
1964!--    forbidden.
1965       IF ( ibc_sa_t == 0  .AND.   constant_top_salinityflux  .AND. &
1966            top_salinityflux /= 0.0_wp )  THEN
1967          message_string = 'boundary condition: bc_sa_t = "' // &
1968                           TRIM( bc_sa_t ) // '" is not allowed with ' // &
1969                           'constant_top_salinityflux = .TRUE.'
1970          CALL message( 'check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
1971       ENDIF
1972
1973    ENDIF
1974
1975!
1976!-- In case of humidity or passive scalar, set boundary conditions for total
1977!-- water content / scalar
1978    IF ( humidity  .OR.  passive_scalar ) THEN
1979       IF ( humidity )  THEN
1980          sq = 'q'
1981       ELSE
1982          sq = 's'
1983       ENDIF
1984       IF ( bc_q_b == 'dirichlet' )  THEN
1985          ibc_q_b = 0
1986       ELSEIF ( bc_q_b == 'neumann' )  THEN
1987          ibc_q_b = 1
1988       ELSE
1989          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &
1990                           '_b ="' // TRIM( bc_q_b ) // '"'
1991          CALL message( 'check_parameters', 'PA0071', 1, 2, 0, 6, 0 )
1992       ENDIF
1993       IF ( bc_q_t == 'dirichlet' )  THEN
1994          ibc_q_t = 0
1995       ELSEIF ( bc_q_t == 'neumann' )  THEN
1996          ibc_q_t = 1
1997       ELSE
1998          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &
1999                           '_t ="' // TRIM( bc_q_t ) // '"'
2000          CALL message( 'check_parameters', 'PA0072', 1, 2, 0, 6, 0 )
2001       ENDIF
2002
2003       IF ( surface_waterflux == 9999999.9_wp  )  THEN
2004          constant_waterflux = .FALSE.
2005          IF ( large_scale_forcing )  THEN
2006             IF ( ibc_q_b == 0 )  THEN
2007                constant_waterflux = .FALSE.
2008             ELSEIF ( ibc_q_b == 1 )  THEN
2009                constant_waterflux = .TRUE.
2010                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
2011                   surface_waterflux = qsws_surf(1)
2012                ENDIF
2013             ENDIF
2014          ENDIF
2015       ELSE
2016          constant_waterflux = .TRUE.
2017          IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
2018                 large_scale_forcing ) THEN
2019             surface_waterflux = qsws_surf(1)
2020          ENDIF
2021       ENDIF
2022
2023!
2024!--    A given surface humidity implies Dirichlet boundary condition for
2025!--    humidity. In this case specification of a constant water flux is
2026!--    forbidden.
2027       IF ( ibc_q_b == 0  .AND.  constant_waterflux )  THEN
2028          message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b ' // &
2029                           '= "' // TRIM( bc_q_b ) // '" is not allowed wi' // &
2030                           'th prescribed surface flux'
2031          CALL message( 'check_parameters', 'PA0073', 1, 2, 0, 6, 0 )
2032       ENDIF
2033       IF ( constant_waterflux  .AND.  q_surface_initial_change /= 0.0_wp )  THEN
2034          WRITE( message_string, * )  'a prescribed surface flux is not allo', &
2035                 'wed with ', sq, '_surface_initial_change (/=0) = ', &
2036                 q_surface_initial_change
2037          CALL message( 'check_parameters', 'PA0074', 1, 2, 0, 6, 0 )
2038       ENDIF
2039
2040    ENDIF
2041!
2042!-- Boundary conditions for horizontal components of wind speed
2043    IF ( bc_uv_b == 'dirichlet' )  THEN
2044       ibc_uv_b = 0
2045    ELSEIF ( bc_uv_b == 'neumann' )  THEN
2046       ibc_uv_b = 1
2047       IF ( prandtl_layer )  THEN
2048          message_string = 'boundary condition: bc_uv_b = "' // &
2049               TRIM( bc_uv_b ) // '" is not allowed with prandtl_layer = .TRUE.'
2050          CALL message( 'check_parameters', 'PA0075', 1, 2, 0, 6, 0 )
2051       ENDIF
2052    ELSE
2053       message_string = 'unknown boundary condition: bc_uv_b = "' // &
2054                        TRIM( bc_uv_b ) // '"'
2055       CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
2056    ENDIF
2057!
2058!-- In case of coupled simulations u and v at the ground in atmosphere will be
2059!-- assigned with the u and v values of the ocean surface
2060    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
2061       ibc_uv_b = 2
2062    ENDIF
2063
2064    IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
2065       bc_uv_t = 'neumann'
2066       ibc_uv_t = 1
2067    ELSE
2068       IF ( bc_uv_t == 'dirichlet' .OR. bc_uv_t == 'dirichlet_0' )  THEN
2069          ibc_uv_t = 0
2070          IF ( bc_uv_t == 'dirichlet_0' )  THEN
2071!
2072!--          Velocities for the initial u,v-profiles are set zero at the top
2073!--          in case of dirichlet_0 conditions
2074             u_init(nzt+1)    = 0.0_wp
2075             v_init(nzt+1)    = 0.0_wp
2076          ENDIF
2077       ELSEIF ( bc_uv_t == 'neumann' )  THEN
2078          ibc_uv_t = 1
2079       ELSE
2080          message_string = 'unknown boundary condition: bc_uv_t = "' // &
2081                           TRIM( bc_uv_t ) // '"'
2082          CALL message( 'check_parameters', 'PA0077', 1, 2, 0, 6, 0 )
2083       ENDIF
2084    ENDIF
2085
2086!
2087!-- Compute and check, respectively, the Rayleigh Damping parameter
2088    IF ( rayleigh_damping_factor == -1.0_wp )  THEN
2089       rayleigh_damping_factor = 0.0_wp
2090    ELSE
2091       IF ( rayleigh_damping_factor < 0.0_wp .OR. rayleigh_damping_factor > 1.0_wp ) &
2092       THEN
2093          WRITE( message_string, * )  'rayleigh_damping_factor = ', &
2094                              rayleigh_damping_factor, ' out of range [0.0,1.0]'
2095          CALL message( 'check_parameters', 'PA0078', 1, 2, 0, 6, 0 )
2096       ENDIF
2097    ENDIF
2098
2099    IF ( rayleigh_damping_height == -1.0_wp )  THEN
2100       IF ( .NOT. ocean )  THEN
2101          rayleigh_damping_height = 0.66666666666_wp * zu(nzt)
2102       ELSE
2103          rayleigh_damping_height = 0.66666666666_wp * zu(nzb)
2104       ENDIF
2105    ELSE
2106       IF ( .NOT. ocean )  THEN
2107          IF ( rayleigh_damping_height < 0.0_wp  .OR. &
2108               rayleigh_damping_height > zu(nzt) )  THEN
2109             WRITE( message_string, * )  'rayleigh_damping_height = ', &
2110                   rayleigh_damping_height, ' out of range [0.0,', zu(nzt), ']'
2111             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
2112          ENDIF
2113       ELSE
2114          IF ( rayleigh_damping_height > 0.0_wp  .OR. &
2115               rayleigh_damping_height < zu(nzb) )  THEN
2116             WRITE( message_string, * )  'rayleigh_damping_height = ', &
2117                   rayleigh_damping_height, ' out of range [0.0,', zu(nzb), ']'
2118             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
2119          ENDIF
2120       ENDIF
2121    ENDIF
2122
2123!
2124!-- Check number of chosen statistic regions. More than 10 regions are not
2125!-- allowed, because so far no more than 10 corresponding output files can
2126!-- be opened (cf. check_open)
2127    IF ( statistic_regions > 9  .OR.  statistic_regions < 0 )  THEN
2128       WRITE ( message_string, * ) 'number of statistic_regions = ', &
2129                   statistic_regions+1, ' but only 10 regions are allowed'
2130       CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 )
2131    ENDIF
2132    IF ( normalizing_region > statistic_regions  .OR. &
2133         normalizing_region < 0)  THEN
2134       WRITE ( message_string, * ) 'normalizing_region = ', &
2135                normalizing_region, ' must be >= 0 and <= ',statistic_regions, &
2136                ' (value of statistic_regions)'
2137       CALL message( 'check_parameters', 'PA0083', 1, 2, 0, 6, 0 )
2138    ENDIF
2139
2140!
2141!-- Set the default intervals for data output, if necessary
2142!-- NOTE: dt_dosp has already been set in package_parin
2143    IF ( dt_data_output /= 9999999.9_wp )  THEN
2144       IF ( dt_dopr           == 9999999.9_wp )  dt_dopr           = dt_data_output
2145       IF ( dt_dopts          == 9999999.9_wp )  dt_dopts          = dt_data_output
2146       IF ( dt_do2d_xy        == 9999999.9_wp )  dt_do2d_xy        = dt_data_output
2147       IF ( dt_do2d_xz        == 9999999.9_wp )  dt_do2d_xz        = dt_data_output
2148       IF ( dt_do2d_yz        == 9999999.9_wp )  dt_do2d_yz        = dt_data_output
2149       IF ( dt_do3d           == 9999999.9_wp )  dt_do3d           = dt_data_output
2150       IF ( dt_data_output_av == 9999999.9_wp )  dt_data_output_av = dt_data_output
2151       DO  mid = 1, max_masks
2152          IF ( dt_domask(mid) == 9999999.9_wp )  dt_domask(mid)    = dt_data_output
2153       ENDDO
2154    ENDIF
2155
2156!
2157!-- Set the default skip time intervals for data output, if necessary
2158    IF ( skip_time_dopr    == 9999999.9_wp ) &
2159                                       skip_time_dopr    = skip_time_data_output
2160    IF ( skip_time_dosp    == 9999999.9_wp ) &
2161                                       skip_time_dosp    = skip_time_data_output
2162    IF ( skip_time_do2d_xy == 9999999.9_wp ) &
2163                                       skip_time_do2d_xy = skip_time_data_output
2164    IF ( skip_time_do2d_xz == 9999999.9_wp ) &
2165                                       skip_time_do2d_xz = skip_time_data_output
2166    IF ( skip_time_do2d_yz == 9999999.9_wp ) &
2167                                       skip_time_do2d_yz = skip_time_data_output
2168    IF ( skip_time_do3d    == 9999999.9_wp ) &
2169                                       skip_time_do3d    = skip_time_data_output
2170    IF ( skip_time_data_output_av == 9999999.9_wp ) &
2171                                skip_time_data_output_av = skip_time_data_output
2172    DO  mid = 1, max_masks
2173       IF ( skip_time_domask(mid) == 9999999.9_wp ) &
2174                                skip_time_domask(mid)    = skip_time_data_output
2175    ENDDO
2176
2177!
2178!-- Check the average intervals (first for 3d-data, then for profiles and
2179!-- spectra)
2180    IF ( averaging_interval > dt_data_output_av )  THEN
2181       WRITE( message_string, * )  'averaging_interval = ', &
2182             averaging_interval, ' must be <= dt_data_output = ', dt_data_output
2183       CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 )
2184    ENDIF
2185
2186    IF ( averaging_interval_pr == 9999999.9_wp )  THEN
2187       averaging_interval_pr = averaging_interval
2188    ENDIF
2189
2190    IF ( averaging_interval_pr > dt_dopr )  THEN
2191       WRITE( message_string, * )  'averaging_interval_pr = ', &
2192             averaging_interval_pr, ' must be <= dt_dopr = ', dt_dopr
2193       CALL message( 'check_parameters', 'PA0086', 1, 2, 0, 6, 0 )
2194    ENDIF
2195
2196    IF ( averaging_interval_sp == 9999999.9_wp )  THEN
2197       averaging_interval_sp = averaging_interval
2198    ENDIF
2199
2200    IF ( averaging_interval_sp > dt_dosp )  THEN
2201       WRITE( message_string, * )  'averaging_interval_sp = ', &
2202             averaging_interval_sp, ' must be <= dt_dosp = ', dt_dosp
2203       CALL message( 'check_parameters', 'PA0087', 1, 2, 0, 6, 0 )
2204    ENDIF
2205
2206!
2207!-- Set the default interval for profiles entering the temporal average
2208    IF ( dt_averaging_input_pr == 9999999.9_wp )  THEN
2209       dt_averaging_input_pr = dt_averaging_input
2210    ENDIF
2211
2212!
2213!-- Set the default interval for the output of timeseries to a reasonable
2214!-- value (tries to minimize the number of calls of flow_statistics)
2215    IF ( dt_dots == 9999999.9_wp )  THEN
2216       IF ( averaging_interval_pr == 0.0_wp )  THEN
2217          dt_dots = MIN( dt_run_control, dt_dopr )
2218       ELSE
2219          dt_dots = MIN( dt_run_control, dt_averaging_input_pr )
2220       ENDIF
2221    ENDIF
2222
2223!
2224!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
2225    IF ( dt_averaging_input > averaging_interval )  THEN
2226       WRITE( message_string, * )  'dt_averaging_input = ', &
2227                dt_averaging_input, ' must be <= averaging_interval = ', &
2228                averaging_interval
2229       CALL message( 'check_parameters', 'PA0088', 1, 2, 0, 6, 0 )
2230    ENDIF
2231
2232    IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
2233       WRITE( message_string, * )  'dt_averaging_input_pr = ', &
2234                dt_averaging_input_pr, ' must be <= averaging_interval_pr = ', &
2235                averaging_interval_pr
2236       CALL message( 'check_parameters', 'PA0089', 1, 2, 0, 6, 0 )
2237    ENDIF
2238
2239!
2240!-- Set the default value for the integration interval of precipitation amount
2241    IF ( precipitation )  THEN
2242       IF ( precipitation_amount_interval == 9999999.9_wp )  THEN
2243          precipitation_amount_interval = dt_do2d_xy
2244       ELSE
2245          IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
2246             WRITE( message_string, * )  'precipitation_amount_interval = ', &
2247                 precipitation_amount_interval, ' must not be larger than ', &
2248                 'dt_do2d_xy = ', dt_do2d_xy
2249             CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
2250          ENDIF
2251       ENDIF
2252    ENDIF
2253
2254!
2255!-- Determine the number of output profiles and check whether they are
2256!-- permissible
2257    DO  WHILE ( data_output_pr(dopr_n+1) /= '          ' )
2258
2259       dopr_n = dopr_n + 1
2260       i = dopr_n
2261
2262!
2263!--    Determine internal profile number (for hom, homs)
2264!--    and store height levels
2265       SELECT CASE ( TRIM( data_output_pr(i) ) )
2266
2267          CASE ( 'u', '#u' )
2268             dopr_index(i) = 1
2269             dopr_unit(i)  = 'm/s'
2270             hom(:,2,1,:)  = SPREAD( zu, 2, statistic_regions+1 )
2271             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2272                dopr_initial_index(i) = 5
2273                hom(:,2,5,:)          = SPREAD( zu, 2, statistic_regions+1 )
2274                data_output_pr(i)     = data_output_pr(i)(2:)
2275             ENDIF
2276
2277          CASE ( 'v', '#v' )
2278             dopr_index(i) = 2
2279             dopr_unit(i)  = 'm/s'
2280             hom(:,2,2,:)  = SPREAD( zu, 2, statistic_regions+1 )
2281             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2282                dopr_initial_index(i) = 6
2283                hom(:,2,6,:)          = SPREAD( zu, 2, statistic_regions+1 )
2284                data_output_pr(i)     = data_output_pr(i)(2:)
2285             ENDIF
2286
2287          CASE ( 'w' )
2288             dopr_index(i) = 3
2289             dopr_unit(i)  = 'm/s'
2290             hom(:,2,3,:)  = SPREAD( zw, 2, statistic_regions+1 )
2291
2292          CASE ( 'pt', '#pt' )
2293             IF ( .NOT. cloud_physics ) THEN
2294                dopr_index(i) = 4
2295                dopr_unit(i)  = 'K'
2296                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2297                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2298                   dopr_initial_index(i) = 7
2299                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2300                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2301                   data_output_pr(i)     = data_output_pr(i)(2:)
2302                ENDIF
2303             ELSE
2304                dopr_index(i) = 43
2305                dopr_unit(i)  = 'K'
2306                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
2307                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2308                   dopr_initial_index(i) = 28
2309                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
2310                   hom(nzb,2,28,:)       = 0.0_wp    ! because zu(nzb) is negative
2311                   data_output_pr(i)     = data_output_pr(i)(2:)
2312                ENDIF
2313             ENDIF
2314
2315          CASE ( 'e' )
2316             dopr_index(i)  = 8
2317             dopr_unit(i)   = 'm2/s2'
2318             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
2319             hom(nzb,2,8,:) = 0.0_wp
2320
2321          CASE ( 'km', '#km' )
2322             dopr_index(i)  = 9
2323             dopr_unit(i)   = 'm2/s'
2324             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
2325             hom(nzb,2,9,:) = 0.0_wp
2326             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2327                dopr_initial_index(i) = 23
2328                hom(:,2,23,:)         = hom(:,2,9,:)
2329                data_output_pr(i)     = data_output_pr(i)(2:)
2330             ENDIF
2331
2332          CASE ( 'kh', '#kh' )
2333             dopr_index(i)   = 10
2334             dopr_unit(i)    = 'm2/s'
2335             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
2336             hom(nzb,2,10,:) = 0.0_wp
2337             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2338                dopr_initial_index(i) = 24
2339                hom(:,2,24,:)         = hom(:,2,10,:)
2340                data_output_pr(i)     = data_output_pr(i)(2:)
2341             ENDIF
2342
2343          CASE ( 'l', '#l' )
2344             dopr_index(i)   = 11
2345             dopr_unit(i)    = 'm'
2346             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
2347             hom(nzb,2,11,:) = 0.0_wp
2348             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2349                dopr_initial_index(i) = 25
2350                hom(:,2,25,:)         = hom(:,2,11,:)
2351                data_output_pr(i)     = data_output_pr(i)(2:)
2352             ENDIF
2353
2354          CASE ( 'w"u"' )
2355             dopr_index(i) = 12
2356             dopr_unit(i)  = 'm2/s2'
2357             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
2358             IF ( prandtl_layer )  hom(nzb,2,12,:) = zu(1)
2359
2360          CASE ( 'w*u*' )
2361             dopr_index(i) = 13
2362             dopr_unit(i)  = 'm2/s2'
2363             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
2364
2365          CASE ( 'w"v"' )
2366             dopr_index(i) = 14
2367             dopr_unit(i)  = 'm2/s2'
2368             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
2369             IF ( prandtl_layer )  hom(nzb,2,14,:) = zu(1)
2370
2371          CASE ( 'w*v*' )
2372             dopr_index(i) = 15
2373             dopr_unit(i)  = 'm2/s2'
2374             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
2375
2376          CASE ( 'w"pt"' )
2377             dopr_index(i) = 16
2378             dopr_unit(i)  = 'K m/s'
2379             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
2380
2381          CASE ( 'w*pt*' )
2382             dopr_index(i) = 17
2383             dopr_unit(i)  = 'K m/s'
2384             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
2385
2386          CASE ( 'wpt' )
2387             dopr_index(i) = 18
2388             dopr_unit(i)  = 'K m/s'
2389             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
2390
2391          CASE ( 'wu' )
2392             dopr_index(i) = 19
2393             dopr_unit(i)  = 'm2/s2'
2394             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
2395             IF ( prandtl_layer )  hom(nzb,2,19,:) = zu(1)
2396
2397          CASE ( 'wv' )
2398             dopr_index(i) = 20
2399             dopr_unit(i)  = 'm2/s2'
2400             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
2401             IF ( prandtl_layer )  hom(nzb,2,20,:) = zu(1)
2402
2403          CASE ( 'w*pt*BC' )
2404             dopr_index(i) = 21
2405             dopr_unit(i)  = 'K m/s'
2406             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
2407
2408          CASE ( 'wptBC' )
2409             dopr_index(i) = 22
2410             dopr_unit(i)  = 'K m/s'
2411             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
2412
2413          CASE ( 'sa', '#sa' )
2414             IF ( .NOT. ocean )  THEN
2415                message_string = 'data_output_pr = ' // &
2416                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2417                                 'lemented for ocean = .FALSE.'
2418                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2419             ELSE
2420                dopr_index(i) = 23
2421                dopr_unit(i)  = 'psu'
2422                hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 )
2423                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2424                   dopr_initial_index(i) = 26
2425                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2426                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2427                   data_output_pr(i)     = data_output_pr(i)(2:)
2428                ENDIF
2429             ENDIF
2430
2431          CASE ( 'u*2' )
2432             dopr_index(i) = 30
2433             dopr_unit(i)  = 'm2/s2'
2434             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
2435
2436          CASE ( 'v*2' )
2437             dopr_index(i) = 31
2438             dopr_unit(i)  = 'm2/s2'
2439             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
2440
2441          CASE ( 'w*2' )
2442             dopr_index(i) = 32
2443             dopr_unit(i)  = 'm2/s2'
2444             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
2445
2446          CASE ( 'pt*2' )
2447             dopr_index(i) = 33
2448             dopr_unit(i)  = 'K2'
2449             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
2450
2451          CASE ( 'e*' )
2452             dopr_index(i) = 34
2453             dopr_unit(i)  = 'm2/s2'
2454             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
2455
2456          CASE ( 'w*2pt*' )
2457             dopr_index(i) = 35
2458             dopr_unit(i)  = 'K m2/s2'
2459             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
2460
2461          CASE ( 'w*pt*2' )
2462             dopr_index(i) = 36
2463             dopr_unit(i)  = 'K2 m/s'
2464             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
2465
2466          CASE ( 'w*e*' )
2467             dopr_index(i) = 37
2468             dopr_unit(i)  = 'm3/s3'
2469             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
2470
2471          CASE ( 'w*3' )
2472             dopr_index(i) = 38
2473             dopr_unit(i)  = 'm3/s3'
2474             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
2475
2476          CASE ( 'Sw' )
2477             dopr_index(i) = 39
2478             dopr_unit(i)  = 'none'
2479             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
2480
2481          CASE ( 'p' )
2482             dopr_index(i) = 40
2483             dopr_unit(i)  = 'Pa'
2484             hom(:,2,40,:) = SPREAD( zu, 2, statistic_regions+1 )
2485
2486          CASE ( 'q', '#q' )
2487             IF ( .NOT. humidity )  THEN
2488                message_string = 'data_output_pr = ' // &
2489                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2490                                 'lemented for humidity = .FALSE.'
2491                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2492             ELSE
2493                dopr_index(i) = 41
2494                dopr_unit(i)  = 'kg/kg'
2495                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2496                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2497                   dopr_initial_index(i) = 26
2498                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2499                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2500                   data_output_pr(i)     = data_output_pr(i)(2:)
2501                ENDIF
2502             ENDIF
2503
2504          CASE ( 's', '#s' )
2505             IF ( .NOT. passive_scalar )  THEN
2506                message_string = 'data_output_pr = ' // &
2507                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2508                                 'lemented for passive_scalar = .FALSE.'
2509                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2510             ELSE
2511                dopr_index(i) = 41
2512                dopr_unit(i)  = 'kg/m3'
2513                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2514                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2515                   dopr_initial_index(i) = 26
2516                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2517                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2518                   data_output_pr(i)     = data_output_pr(i)(2:)
2519                ENDIF
2520             ENDIF
2521
2522          CASE ( 'qv', '#qv' )
2523             IF ( .NOT. cloud_physics ) THEN
2524                dopr_index(i) = 41
2525                dopr_unit(i)  = 'kg/kg'
2526                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2527                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2528                   dopr_initial_index(i) = 26
2529                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2530                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2531                   data_output_pr(i)     = data_output_pr(i)(2:)
2532                ENDIF
2533             ELSE
2534                dopr_index(i) = 42
2535                dopr_unit(i)  = 'kg/kg'
2536                hom(:,2,42,:) = SPREAD( zu, 2, statistic_regions+1 )
2537                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2538                   dopr_initial_index(i) = 27
2539                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
2540                   hom(nzb,2,27,:)       = 0.0_wp   ! because zu(nzb) is negative
2541                   data_output_pr(i)     = data_output_pr(i)(2:)
2542                ENDIF
2543             ENDIF
2544
2545          CASE ( 'lpt', '#lpt' )
2546             IF ( .NOT. cloud_physics ) THEN
2547                message_string = 'data_output_pr = ' // &
2548                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2549                                 'lemented for cloud_physics = .FALSE.'
2550                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2551             ELSE
2552                dopr_index(i) = 4
2553                dopr_unit(i)  = 'K'
2554                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2555                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2556                   dopr_initial_index(i) = 7
2557                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2558                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2559                   data_output_pr(i)     = data_output_pr(i)(2:)
2560                ENDIF
2561             ENDIF
2562
2563          CASE ( 'vpt', '#vpt' )
2564             dopr_index(i) = 44
2565             dopr_unit(i)  = 'K'
2566             hom(:,2,44,:) = SPREAD( zu, 2, statistic_regions+1 )
2567             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2568                dopr_initial_index(i) = 29
2569                hom(:,2,29,:)         = SPREAD( zu, 2, statistic_regions+1 )
2570                hom(nzb,2,29,:)       = 0.0_wp    ! because zu(nzb) is negative
2571                data_output_pr(i)     = data_output_pr(i)(2:)
2572             ENDIF
2573
2574          CASE ( 'w"vpt"' )
2575             dopr_index(i) = 45
2576             dopr_unit(i)  = 'K m/s'
2577             hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
2578
2579          CASE ( 'w*vpt*' )
2580             dopr_index(i) = 46
2581             dopr_unit(i)  = 'K m/s'
2582             hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
2583
2584          CASE ( 'wvpt' )
2585             dopr_index(i) = 47
2586             dopr_unit(i)  = 'K m/s'
2587             hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
2588
2589          CASE ( 'w"q"' )
2590             IF ( .NOT. humidity )  THEN
2591                message_string = 'data_output_pr = ' // &
2592                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2593                                 'lemented for humidity = .FALSE.'
2594                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2595             ELSE
2596                dopr_index(i) = 48
2597                dopr_unit(i)  = 'kg/kg m/s'
2598                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2599             ENDIF
2600
2601          CASE ( 'w*q*' )
2602             IF ( .NOT. humidity )  THEN
2603                message_string = 'data_output_pr = ' // &
2604                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2605                                 'lemented for humidity = .FALSE.'
2606                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2607             ELSE
2608                dopr_index(i) = 49
2609                dopr_unit(i)  = 'kg/kg m/s'
2610                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2611             ENDIF
2612
2613          CASE ( 'wq' )
2614             IF ( .NOT. humidity )  THEN
2615                message_string = 'data_output_pr = ' // &
2616                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2617                                 'lemented for humidity = .FALSE.'
2618                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2619             ELSE
2620                dopr_index(i) = 50
2621                dopr_unit(i)  = 'kg/kg m/s'
2622                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2623             ENDIF
2624
2625          CASE ( 'w"s"' )
2626             IF ( .NOT. passive_scalar ) THEN
2627                message_string = 'data_output_pr = ' // &
2628                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2629                                 'lemented for passive_scalar = .FALSE.'
2630                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2631             ELSE
2632                dopr_index(i) = 48
2633                dopr_unit(i)  = 'kg/m3 m/s'
2634                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2635             ENDIF
2636
2637          CASE ( 'w*s*' )
2638             IF ( .NOT. passive_scalar ) THEN
2639                message_string = 'data_output_pr = ' // &
2640                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2641                                 'lemented for passive_scalar = .FALSE.'
2642                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2643             ELSE
2644                dopr_index(i) = 49
2645                dopr_unit(i)  = 'kg/m3 m/s'
2646                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2647             ENDIF
2648
2649          CASE ( 'ws' )
2650             IF ( .NOT. passive_scalar ) THEN
2651                message_string = 'data_output_pr = ' // &
2652                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2653                                 'lemented for passive_scalar = .FALSE.'
2654                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2655             ELSE
2656                dopr_index(i) = 50
2657                dopr_unit(i)  = 'kg/m3 m/s'
2658                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2659             ENDIF
2660
2661          CASE ( 'w"qv"' )
2662             IF ( humidity  .AND.  .NOT. cloud_physics ) &
2663             THEN
2664                dopr_index(i) = 48
2665                dopr_unit(i)  = 'kg/kg m/s'
2666                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2667             ELSEIF( humidity .AND. cloud_physics ) THEN
2668                dopr_index(i) = 51
2669                dopr_unit(i)  = 'kg/kg m/s'
2670                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
2671             ELSE
2672                message_string = 'data_output_pr = ' // &
2673                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2674                                 'lemented for cloud_physics = .FALSE. an&' // &
2675                                 'd humidity = .FALSE.'
2676                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2677             ENDIF
2678
2679          CASE ( 'w*qv*' )
2680             IF ( humidity  .AND.  .NOT. cloud_physics ) &
2681             THEN
2682                dopr_index(i) = 49
2683                dopr_unit(i)  = 'kg/kg m/s'
2684                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2685             ELSEIF( humidity .AND. cloud_physics ) THEN
2686                dopr_index(i) = 52
2687                dopr_unit(i)  = 'kg/kg m/s'
2688                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
2689             ELSE
2690                message_string = 'data_output_pr = ' // &
2691                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2692                                 'lemented for cloud_physics = .FALSE. an&' // &
2693                                 'd humidity = .FALSE.'
2694                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2695             ENDIF
2696
2697          CASE ( 'wqv' )
2698             IF ( humidity  .AND.  .NOT. cloud_physics ) &
2699             THEN
2700                dopr_index(i) = 50
2701                dopr_unit(i)  = 'kg/kg m/s'
2702                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2703             ELSEIF( humidity .AND. cloud_physics ) THEN
2704                dopr_index(i) = 53
2705                dopr_unit(i)  = 'kg/kg m/s'
2706                hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
2707             ELSE
2708                message_string = 'data_output_pr = ' //                        &
2709                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2710                                 'lemented for cloud_physics = .FALSE. an&' // &
2711                                 'd humidity = .FALSE.'
2712                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2713             ENDIF
2714
2715          CASE ( 'ql' )
2716             IF ( .NOT. cloud_physics  .AND.  .NOT. cloud_droplets )  THEN
2717                message_string = 'data_output_pr = ' //                        &
2718                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2719                                 'lemented for cloud_physics = .FALSE. or'  // &
2720                                 '&cloud_droplets = .FALSE.'
2721                CALL message( 'check_parameters', 'PA0096', 1, 2, 0, 6, 0 )
2722             ELSE
2723                dopr_index(i) = 54
2724                dopr_unit(i)  = 'kg/kg'
2725                hom(:,2,54,:)  = SPREAD( zu, 2, statistic_regions+1 )
2726             ENDIF
2727
2728          CASE ( 'w*u*u*:dz' )
2729             dopr_index(i) = 55
2730             dopr_unit(i)  = 'm2/s3'
2731             hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
2732
2733          CASE ( 'w*p*:dz' )
2734             dopr_index(i) = 56
2735             dopr_unit(i)  = 'm2/s3'
2736             hom(:,2,56,:) = SPREAD( zw, 2, statistic_regions+1 )
2737
2738          CASE ( 'w"e:dz' )
2739             dopr_index(i) = 57
2740             dopr_unit(i)  = 'm2/s3'
2741             hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
2742
2743
2744          CASE ( 'u"pt"' )
2745             dopr_index(i) = 58
2746             dopr_unit(i)  = 'K m/s'
2747             hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
2748
2749          CASE ( 'u*pt*' )
2750             dopr_index(i) = 59
2751             dopr_unit(i)  = 'K m/s'
2752             hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
2753
2754          CASE ( 'upt_t' )
2755             dopr_index(i) = 60
2756             dopr_unit(i)  = 'K m/s'
2757             hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
2758
2759          CASE ( 'v"pt"' )
2760             dopr_index(i) = 61
2761             dopr_unit(i)  = 'K m/s'
2762             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
2763             
2764          CASE ( 'v*pt*' )
2765             dopr_index(i) = 62
2766             dopr_unit(i)  = 'K m/s'
2767             hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
2768
2769          CASE ( 'vpt_t' )
2770             dopr_index(i) = 63
2771             dopr_unit(i)  = 'K m/s'
2772             hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
2773
2774          CASE ( 'rho' )
2775             IF ( .NOT. ocean ) THEN
2776                message_string = 'data_output_pr = ' //                        &
2777                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2778                                 'lemented for ocean = .FALSE.'
2779                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2780             ELSE
2781                dopr_index(i) = 64
2782                dopr_unit(i)  = 'kg/m3'
2783                hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
2784                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2785                   dopr_initial_index(i) = 77
2786                   hom(:,2,77,:)         = SPREAD( zu, 2, statistic_regions+1 )
2787                   hom(nzb,2,77,:)       = 0.0_wp    ! because zu(nzb) is negative
2788                   data_output_pr(i)     = data_output_pr(i)(2:)
2789                ENDIF
2790             ENDIF
2791
2792          CASE ( 'w"sa"' )
2793             IF ( .NOT. ocean ) THEN
2794                message_string = 'data_output_pr = ' //                        &
2795                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2796                                 'lemented for ocean = .FALSE.'
2797                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2798             ELSE
2799                dopr_index(i) = 65
2800                dopr_unit(i)  = 'psu m/s'
2801                hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 )
2802             ENDIF
2803
2804          CASE ( 'w*sa*' )
2805             IF ( .NOT. ocean ) THEN
2806                message_string = 'data_output_pr = ' //                        &
2807                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2808                                 'lemented for ocean = .FALSE.'
2809                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2810             ELSE
2811                dopr_index(i) = 66
2812                dopr_unit(i)  = 'psu m/s'
2813                hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 )
2814             ENDIF
2815
2816          CASE ( 'wsa' )
2817             IF ( .NOT. ocean ) THEN
2818                message_string = 'data_output_pr = ' //                        &
2819                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2820                                 'lemented for ocean = .FALSE.'
2821                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2822             ELSE
2823                dopr_index(i) = 67
2824                dopr_unit(i)  = 'psu m/s'
2825                hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
2826             ENDIF
2827
2828          CASE ( 'w*p*' )
2829             dopr_index(i) = 68
2830             dopr_unit(i)  = 'm3/s3'
2831             hom(:,2,68,:) = SPREAD( zu, 2, statistic_regions+1 )
2832
2833          CASE ( 'w"e' )
2834             dopr_index(i) = 69
2835             dopr_unit(i)  = 'm3/s3'
2836             hom(:,2,69,:) = SPREAD( zu, 2, statistic_regions+1 )
2837
2838          CASE ( 'q*2' )
2839             IF ( .NOT. humidity )  THEN
2840                message_string = 'data_output_pr = ' //                        &
2841                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2842                                 'lemented for humidity = .FALSE.'
2843                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2844             ELSE
2845                dopr_index(i) = 70
2846                dopr_unit(i)  = 'kg2/kg2'
2847                hom(:,2,70,:) = SPREAD( zu, 2, statistic_regions+1 )
2848             ENDIF
2849
2850          CASE ( 'prho' )
2851             IF ( .NOT. ocean ) THEN
2852                message_string = 'data_output_pr = ' //                        &
2853                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2854                                 'lemented for ocean = .FALSE.'
2855                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2856             ELSE
2857                dopr_index(i) = 71
2858                dopr_unit(i)  = 'kg/m3'
2859                hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
2860             ENDIF
2861
2862          CASE ( 'hyp' )
2863             dopr_index(i) = 72
2864             dopr_unit(i)  = 'dbar'
2865             hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
2866
2867          CASE ( 'nr' )
2868             IF ( .NOT. cloud_physics )  THEN
2869                message_string = 'data_output_pr = ' //                        &
2870                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2871                                 'lemented for cloud_physics = .FALSE.'
2872                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2873             ELSEIF ( icloud_scheme /= 0 )  THEN
2874                message_string = 'data_output_pr = ' //                        &
2875                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2876                                 'lemented for cloud_scheme /= seifert_beheng'
2877                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2878             ELSEIF ( .NOT. precipitation )  THEN
2879                message_string = 'data_output_pr = ' //                        &
2880                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2881                                 'lemented for precipitation = .FALSE.'
2882                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
2883             ELSE
2884                dopr_index(i) = 73
2885                dopr_unit(i)  = '1/m3'
2886                hom(:,2,73,:)  = SPREAD( zu, 2, statistic_regions+1 )
2887             ENDIF
2888
2889          CASE ( 'qr' )
2890             IF ( .NOT. cloud_physics )  THEN
2891                message_string = 'data_output_pr = ' //                        &
2892                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2893                                 'lemented for cloud_physics = .FALSE.'
2894                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2895             ELSEIF ( icloud_scheme /= 0 )  THEN
2896                message_string = 'data_output_pr = ' //                        &
2897                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2898                                 'lemented for cloud_scheme /= seifert_beheng'
2899                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2900             ELSEIF ( .NOT. precipitation )  THEN
2901                message_string = 'data_output_pr = ' //                        &
2902                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2903                                 'lemented for precipitation = .FALSE.'
2904                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
2905             ELSE
2906                dopr_index(i) = 74
2907                dopr_unit(i)  = 'kg/kg'
2908                hom(:,2,74,:)  = SPREAD( zu, 2, statistic_regions+1 )
2909             ENDIF
2910
2911          CASE ( 'qc' )
2912             IF ( .NOT. cloud_physics )  THEN
2913                message_string = 'data_output_pr = ' //                        &
2914                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2915                                 'lemented for cloud_physics = .FALSE.'
2916                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2917             ELSEIF ( icloud_scheme /= 0 )  THEN
2918                message_string = 'data_output_pr = ' //                        &
2919                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2920                                 'lemented for cloud_scheme /= seifert_beheng'
2921                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2922             ELSE
2923                dopr_index(i) = 75
2924                dopr_unit(i)  = 'kg/kg'
2925                hom(:,2,75,:)  = SPREAD( zu, 2, statistic_regions+1 )
2926             ENDIF
2927
2928          CASE ( 'prr' )
2929             IF ( .NOT. cloud_physics )  THEN
2930                message_string = 'data_output_pr = ' //                        &
2931                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2932                                 'lemented for cloud_physics = .FALSE.'
2933                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2934             ELSEIF ( icloud_scheme /= 0 )  THEN
2935                message_string = 'data_output_pr = ' //                        &
2936                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2937                                 'lemented for cloud_scheme /= seifert_beheng'
2938                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2939             ELSEIF ( .NOT. precipitation )  THEN
2940                message_string = 'data_output_pr = ' //                        &
2941                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2942                                 'lemented for precipitation = .FALSE.'
2943                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
2944
2945             ELSE
2946                dopr_index(i) = 76
2947                dopr_unit(i)  = 'kg/kg m/s'
2948                hom(:,2,76,:)  = SPREAD( zu, 2, statistic_regions+1 )
2949             ENDIF
2950
2951          CASE ( 'ug' )
2952             dopr_index(i) = 78
2953             dopr_unit(i)  = 'm/s'
2954             hom(:,2,78,:) = SPREAD( zu, 2, statistic_regions+1 )
2955
2956          CASE ( 'vg' )
2957             dopr_index(i) = 79
2958             dopr_unit(i)  = 'm/s'
2959             hom(:,2,79,:) = SPREAD( zu, 2, statistic_regions+1 )
2960
2961          CASE ( 'w_subs' )
2962             IF ( .NOT. large_scale_subsidence )  THEN
2963                message_string = 'data_output_pr = ' //                        &
2964                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2965                                 'lemented for large_scale_subsidence = .FALSE.'
2966                CALL message( 'check_parameters', 'PA0382', 1, 2, 0, 6, 0 )
2967             ELSE
2968                dopr_index(i) = 80
2969                dopr_unit(i)  = 'm/s'
2970                hom(:,2,80,:) = SPREAD( zu, 2, statistic_regions+1 )
2971             ENDIF
2972
2973          CASE ( 'td_lsa_lpt' )
2974             IF ( .NOT. large_scale_forcing )  THEN
2975                message_string = 'data_output_pr = ' //                        &
2976                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2977                                 'lemented for large_scale_forcing = .FALSE.'
2978                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
2979             ELSE
2980                dopr_index(i) = 81
2981                dopr_unit(i)  = 'K/s'
2982                hom(:,2,81,:) = SPREAD( zu, 2, statistic_regions+1 )
2983             ENDIF
2984
2985          CASE ( 'td_lsa_q' )
2986             IF ( .NOT. large_scale_forcing )  THEN
2987                message_string = 'data_output_pr = ' //                        &
2988                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2989                                 'lemented for large_scale_forcing = .FALSE.'
2990                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
2991             ELSE
2992                dopr_index(i) = 82
2993                dopr_unit(i)  = 'kg/kgs'
2994                hom(:,2,82,:) = SPREAD( zu, 2, statistic_regions+1 )
2995             ENDIF
2996
2997          CASE ( 'td_sub_lpt' )
2998             IF ( .NOT. large_scale_forcing )  THEN
2999                message_string = 'data_output_pr = ' //                        &
3000                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3001                                 'lemented for large_scale_forcing = .FALSE.'
3002                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
3003             ELSE
3004                dopr_index(i) = 83
3005                dopr_unit(i)  = 'K/s'
3006                hom(:,2,83,:) = SPREAD( zu, 2, statistic_regions+1 )
3007             ENDIF
3008
3009          CASE ( 'td_sub_q' )
3010             IF ( .NOT. large_scale_forcing )  THEN
3011                message_string = 'data_output_pr = ' //                        &
3012                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3013                                 'lemented for large_scale_forcing = .FALSE.'
3014                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
3015             ELSE
3016                dopr_index(i) = 84
3017                dopr_unit(i)  = 'kg/kgs'
3018                hom(:,2,84,:) = SPREAD( zu, 2, statistic_regions+1 )
3019             ENDIF
3020
3021          CASE ( 'td_nud_lpt' )
3022             IF ( .NOT. nudging )  THEN
3023                message_string = 'data_output_pr = ' //                        &
3024                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3025                                 'lemented for nudging = .FALSE.'
3026                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
3027             ELSE
3028                dopr_index(i) = 85
3029                dopr_unit(i)  = 'K/s'
3030                hom(:,2,85,:) = SPREAD( zu, 2, statistic_regions+1 )
3031             ENDIF
3032
3033          CASE ( 'td_nud_q' )
3034             IF ( .NOT. nudging )  THEN
3035                message_string = 'data_output_pr = ' //                        &
3036                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3037                                 'lemented for nudging = .FALSE.'
3038                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
3039             ELSE
3040                dopr_index(i) = 86
3041                dopr_unit(i)  = 'kg/kgs'
3042                hom(:,2,86,:) = SPREAD( zu, 2, statistic_regions+1 )
3043             ENDIF
3044
3045          CASE ( 'td_nud_u' )
3046             IF ( .NOT. nudging )  THEN
3047                message_string = 'data_output_pr = ' //                        &
3048                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3049                                 'lemented for nudging = .FALSE.'
3050                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
3051             ELSE
3052                dopr_index(i) = 87
3053                dopr_unit(i)  = 'm/s2'
3054                hom(:,2,87,:) = SPREAD( zu, 2, statistic_regions+1 )
3055             ENDIF
3056
3057          CASE ( 'td_nud_v' )
3058             IF ( .NOT. nudging )  THEN
3059                message_string = 'data_output_pr = ' //                        &
3060                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3061                                 'lemented for nudging = .FALSE.'
3062                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
3063             ELSE
3064                dopr_index(i) = 88
3065                dopr_unit(i)  = 'm/s2'
3066                hom(:,2,88,:) = SPREAD( zu, 2, statistic_regions+1 )
3067             ENDIF
3068
3069          CASE ( 't_soil', '#t_soil' )
3070             IF ( .NOT. land_surface )  THEN
3071                message_string = 'data_output_pr = ' //                        &
3072                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3073                                 'lemented for land_surface = .FALSE.'
3074                CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
3075             ELSE
3076                dopr_index(i) = 89
3077                dopr_unit(i)  = 'K'
3078                hom(0:nzs-1,2,89,:)  = SPREAD( - zs, 2, statistic_regions+1 )
3079                IF ( data_output_pr(i)(1:1) == '#' )  THEN
3080                   dopr_initial_index(i) = 90
3081                   hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
3082                   data_output_pr(i)     = data_output_pr(i)(2:)
3083                ENDIF
3084             ENDIF
3085
3086          CASE ( 'm_soil', '#m_soil' )
3087             IF ( .NOT. land_surface )  THEN
3088                message_string = 'data_output_pr = ' //                        &
3089                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3090                                 'lemented for land_surface = .FALSE.'
3091                CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
3092             ELSE
3093                dopr_index(i) = 91
3094                dopr_unit(i)  = 'm3/m3'
3095                hom(0:nzs-1,2,91,:)  = SPREAD( - zs, 2, statistic_regions+1 )
3096                IF ( data_output_pr(i)(1:1) == '#' )  THEN
3097                   dopr_initial_index(i) = 92
3098                   hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
3099                   data_output_pr(i)     = data_output_pr(i)(2:)
3100                ENDIF
3101             ENDIF
3102
3103          CASE ( 'rad_lw_in' )
3104             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
3105                message_string = 'data_output_pr = ' //                        &
3106                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3107                                 'lable for radiation = .FALSE. or ' //        &
3108                                 'radiation_scheme = "constant"'
3109                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
3110             ELSE
3111                dopr_index(i) = 102
3112                dopr_unit(i)  = 'W/m2'
3113                hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
3114             ENDIF
3115
3116          CASE ( 'rad_lw_out' )
3117             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
3118                message_string = 'data_output_pr = ' //                        &
3119                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3120                                 'lable for radiation = .FALSE. or ' //        &
3121                                 'radiation_scheme = "constant"'
3122                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
3123             ELSE
3124                dopr_index(i) = 103
3125                dopr_unit(i)  = 'W/m2'
3126                hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
3127             ENDIF
3128
3129          CASE ( 'rad_sw_in' )
3130             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
3131                message_string = 'data_output_pr = ' //                        &
3132                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3133                                 'lable for radiation = .FALSE. or ' //        &
3134                                 'radiation_scheme = "constant"'
3135                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
3136             ELSE
3137                dopr_index(i) = 104
3138                dopr_unit(i)  = 'W/m2'
3139                hom(:,2,104,:)  = SPREAD( zw, 2, statistic_regions+1 )
3140             ENDIF
3141
3142          CASE ( 'rad_sw_out')
3143             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
3144                message_string = 'data_output_pr = ' //                        &
3145                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3146                                 'lable for radiation = .FALSE. or ' //        &
3147                                 'radiation_scheme = "constant"'
3148                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
3149             ELSE
3150                dopr_index(i) = 105
3151                dopr_unit(i)  = 'W/m2'
3152                hom(:,2,105,:)  = SPREAD( zw, 2, statistic_regions+1 )
3153             ENDIF
3154
3155          CASE DEFAULT
3156
3157             CALL user_check_data_output_pr( data_output_pr(i), i, unit )
3158
3159             IF ( unit == 'illegal' )  THEN
3160                IF ( data_output_pr_user(1) /= ' ' )  THEN
3161                   message_string = 'illegal value for data_output_pr or ' //  &
3162                                    'data_output_pr_user = "' //               &
3163                                    TRIM( data_output_pr(i) ) // '"'
3164                   CALL message( 'check_parameters', 'PA0097', 1, 2, 0, 6, 0 )
3165                ELSE
3166                   message_string = 'illegal value for data_output_pr = "' //  &
3167                                    TRIM( data_output_pr(i) ) // '"'
3168                   CALL message( 'check_parameters', 'PA0098', 1, 2, 0, 6, 0 )
3169                ENDIF
3170             ENDIF
3171
3172       END SELECT
3173
3174    ENDDO
3175
3176
3177!
3178!-- Append user-defined data output variables to the standard data output
3179    IF ( data_output_user(1) /= ' ' )  THEN
3180       i = 1
3181       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
3182          i = i + 1
3183       ENDDO
3184       j = 1
3185       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 100 )
3186          IF ( i > 100 )  THEN
3187             message_string = 'number of output quantitities given by data' // &
3188                '_output and data_output_user exceeds the limit of 100'
3189             CALL message( 'check_parameters', 'PA0102', 1, 2, 0, 6, 0 )
3190          ENDIF
3191          data_output(i) = data_output_user(j)
3192          i = i + 1
3193          j = j + 1
3194       ENDDO
3195    ENDIF
3196
3197!
3198!-- Check and set steering parameters for 2d/3d data output and averaging
3199    i   = 1
3200    DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
3201!
3202!--    Check for data averaging
3203       ilen = LEN_TRIM( data_output(i) )
3204       j = 0                                                 ! no data averaging
3205       IF ( ilen > 3 )  THEN
3206          IF ( data_output(i)(ilen-2:ilen) == '_av' )  THEN
3207             j = 1                                           ! data averaging
3208             data_output(i) = data_output(i)(1:ilen-3)
3209          ENDIF
3210       ENDIF
3211!
3212!--    Check for cross section or volume data
3213       ilen = LEN_TRIM( data_output(i) )
3214       k = 0                                                   ! 3d data
3215       var = data_output(i)(1:ilen)
3216       IF ( ilen > 3 )  THEN
3217          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR.                      &
3218               data_output(i)(ilen-2:ilen) == '_xz'  .OR.                      &
3219               data_output(i)(ilen-2:ilen) == '_yz' )  THEN
3220             k = 1                                             ! 2d data
3221             var = data_output(i)(1:ilen-3)
3222          ENDIF
3223       ENDIF
3224
3225!
3226!--    Check for allowed value and set units
3227       SELECT CASE ( TRIM( var ) )
3228
3229          CASE ( 'e' )
3230             IF ( constant_diffusion )  THEN
3231                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3232                                 'res constant_diffusion = .FALSE.'
3233                CALL message( 'check_parameters', 'PA0103', 1, 2, 0, 6, 0 )
3234             ENDIF
3235             unit = 'm2/s2'
3236
3237          CASE ( 'lpt' )
3238             IF ( .NOT. cloud_physics )  THEN
3239                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3240                         'res cloud_physics = .TRUE.'
3241                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3242             ENDIF
3243             unit = 'K'
3244
3245          CASE ( 'm_soil' )
3246             IF ( .NOT. land_surface )  THEN
3247                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3248                         'land_surface = .TRUE.'
3249                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3250             ENDIF
3251             unit = 'm3/m3'
3252
3253          CASE ( 'nr' )
3254             IF ( .NOT. cloud_physics )  THEN
3255                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3256                         'res cloud_physics = .TRUE.'
3257                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3258             ELSEIF ( icloud_scheme /= 0 )  THEN
3259                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3260                         'res cloud_scheme = seifert_beheng'
3261                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3262             ENDIF
3263             unit = '1/m3'
3264
3265          CASE ( 'pc', 'pr' )
3266             IF ( .NOT. particle_advection )  THEN
3267                message_string = 'output of "' // TRIM( var ) // '" requir' // &
3268                   'es a "particles_par"-NAMELIST in the parameter file (PARIN)'
3269                CALL message( 'check_parameters', 'PA0104', 1, 2, 0, 6, 0 )
3270             ENDIF
3271             IF ( TRIM( var ) == 'pc' )  unit = 'number'
3272             IF ( TRIM( var ) == 'pr' )  unit = 'm'
3273
3274          CASE ( 'prr' )
3275             IF ( .NOT. cloud_physics )  THEN
3276                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3277                         'res cloud_physics = .TRUE.'
3278                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3279             ELSEIF ( icloud_scheme /= 0 )  THEN
3280                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3281                         'res cloud_scheme = seifert_beheng'
3282                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3283             ELSEIF ( .NOT. precipitation )  THEN
3284                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3285                                 'res precipitation = .TRUE.'
3286                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3287             ENDIF
3288             unit = 'kg/kg m/s'
3289
3290          CASE ( 'q', 'vpt' )
3291             IF ( .NOT. humidity )  THEN
3292                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3293                                 'res humidity = .TRUE.'
3294                CALL message( 'check_parameters', 'PA0105', 1, 2, 0, 6, 0 )
3295             ENDIF
3296             IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
3297             IF ( TRIM( var ) == 'vpt' )  unit = 'K'
3298
3299          CASE ( 'qc' )
3300             IF ( .NOT. cloud_physics )  THEN
3301                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3302                         'res cloud_physics = .TRUE.'
3303                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3304             ELSEIF ( icloud_scheme /= 0 ) THEN
3305                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3306                         'res cloud_scheme = seifert_beheng'
3307                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3308             ENDIF
3309             unit = 'kg/kg'
3310
3311          CASE ( 'ql' )
3312             IF ( .NOT. ( cloud_physics  .OR.  cloud_droplets ) )  THEN
3313                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3314                         'res cloud_physics = .TRUE. or cloud_droplets = .TRUE.'
3315                CALL message( 'check_parameters', 'PA0106', 1, 2, 0, 6, 0 )
3316             ENDIF
3317             unit = 'kg/kg'
3318
3319          CASE ( 'ql_c', 'ql_v', 'ql_vp' )
3320             IF ( .NOT. cloud_droplets )  THEN
3321                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3322                                 'res cloud_droplets = .TRUE.'
3323                CALL message( 'check_parameters', 'PA0107', 1, 2, 0, 6, 0 )
3324             ENDIF
3325             IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
3326             IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
3327             IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
3328
3329          CASE ( 'qr' )
3330             IF ( .NOT. cloud_physics )  THEN
3331                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3332                         'res cloud_physics = .TRUE.'
3333                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3334             ELSEIF ( icloud_scheme /= 0 ) THEN
3335                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3336                         'res cloud_scheme = seifert_beheng'
3337                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3338             ELSEIF ( .NOT. precipitation )  THEN
3339                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3340                                 'res precipitation = .TRUE.'
3341                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3342             ENDIF
3343             unit = 'kg/kg'
3344
3345          CASE ( 'qv' )
3346             IF ( .NOT. cloud_physics )  THEN
3347                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3348                                 'res cloud_physics = .TRUE.'
3349                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3350             ENDIF
3351             unit = 'kg/kg'
3352
3353
3354          CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out' )
3355             IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' )  THEN
3356                message_string = '"output of "' // TRIM( var ) // '" requi' //  &
3357                                 'res radiation = .TRUE. and ' //              &
3358                                 'radiation_scheme = "rrtmg"'
3359                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
3360             ENDIF
3361             unit = 'W/m2'
3362
3363          CASE ( 'rho' )
3364             IF ( .NOT. ocean )  THEN
3365                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3366                                 'res ocean = .TRUE.'
3367                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
3368             ENDIF
3369             unit = 'kg/m3'
3370
3371          CASE ( 's' )
3372             IF ( .NOT. passive_scalar )  THEN
3373                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3374                                 'res passive_scalar = .TRUE.'
3375                CALL message( 'check_parameters', 'PA0110', 1, 2, 0, 6, 0 )
3376             ENDIF
3377             unit = 'conc'
3378
3379          CASE ( 'sa' )
3380             IF ( .NOT. ocean )  THEN
3381                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3382                                 'res ocean = .TRUE.'
3383                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
3384             ENDIF
3385             unit = 'psu'
3386
3387          CASE ( 't_soil' )
3388             IF ( .NOT. land_surface )  THEN
3389                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3390                         'land_surface = .TRUE.'
3391                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3392             ENDIF
3393             unit = 'K'
3394
3395
3396          CASE ( 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'lai*', 'lwp*',     &
3397                 'm_liq_eb*', 'pra*', 'prr*', 'qsws*', 'qsws_eb*',             &
3398                 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*', 'rad_net*',  &
3399                 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*',   &
3400                 'r_a*', 'r_s*', 'shf*', 'shf_eb*', 't*', 'u*', 'z0*', 'z0h*' )
3401             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
3402                message_string = 'illegal value for data_output: "' //         &
3403                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
3404                                 'cross sections are allowed for this value'
3405                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
3406             ENDIF
3407             IF ( .NOT. radiation .OR. radiation_scheme /= "rrtmg" )  THEN
3408                IF ( TRIM( var ) == 'rrtm_aldif*' .OR.                         &
3409                     TRIM( var ) == 'rrtm_aldir*' .OR.                         &
3410                     TRIM( var ) == 'rrtm_asdif*' .OR.                         &
3411                     TRIM( var ) == 'rrtm_asdir*'      )                       &
3412                THEN
3413                   message_string = 'output of "' // TRIM( var ) // '" require'&
3414                                    // 's radiation = .TRUE. and radiation_sch'&
3415                                    // 'eme = "rrtmg"'
3416                   CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
3417                ENDIF
3418             ENDIF
3419
3420             IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT. land_surface )  THEN
3421                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3422                                 'res land_surface = .TRUE.'
3423                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3424             ENDIF
3425             IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT. land_surface )  THEN
3426                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3427                                 'res land_surface = .TRUE.'
3428                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3429             ENDIF
3430             IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
3431                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3432                                 'res land_surface = .TRUE.'
3433                CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
3434             ENDIF
3435             IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT. land_surface )  THEN
3436                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3437                                 'res land_surface = .TRUE.'
3438                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3439             ENDIF
3440             IF ( TRIM( var ) == 'lai*'  .AND.  .NOT. land_surface )  THEN
3441                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3442                                 'res land_surface = .TRUE.'
3443                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3444             ENDIF
3445             IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
3446                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3447                                 'res cloud_physics = .TRUE.'
3448                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3449             ENDIF
3450             IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT. land_surface )  THEN
3451                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3452                                 'res land_surface = .TRUE.'
3453                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3454             ENDIF
3455             IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
3456                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3457                                 'res precipitation = .TRUE.'
3458                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3459             ENDIF
3460             IF ( TRIM( var ) == 'pra*'  .AND.  j == 1 )  THEN
3461                message_string = 'temporal averaging of precipitation ' //     &
3462                          'amount "' // TRIM( var ) // '" is not possible'
3463                CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
3464             ENDIF
3465             IF ( TRIM( var ) == 'prr*'  .AND.  .NOT. precipitation )  THEN
3466                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3467                                 'res precipitation = .TRUE.'
3468                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3469             ENDIF
3470             IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT. humidity )  THEN
3471                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3472                                 'res humidity = .TRUE.'
3473                CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
3474             ENDIF
3475             IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT. land_surface )  THEN
3476                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3477                                 'res land_surface = .TRUE.'
3478                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3479             ENDIF
3480             IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )  &
3481             THEN
3482                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3483                                 'res land_surface = .TRUE.'
3484                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3485             ENDIF
3486             IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT. land_surface ) &
3487             THEN
3488                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3489                                 'res land_surface = .TRUE.'
3490                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3491             ENDIF
3492             IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )  &
3493             THEN
3494                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3495                                 'res land_surface = .TRUE.'
3496                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3497             ENDIF
3498             IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT. land_surface ) &
3499             THEN
3500                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3501                                 'res land_surface = .TRUE.'
3502                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3503             ENDIF
3504             IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT. land_surface ) &
3505             THEN
3506                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3507                                 'res land_surface = .TRUE.'
3508                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3509             ENDIF
3510
3511             IF ( TRIM( var ) == 'c_liq*' )  unit = 'none'
3512             IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
3513             IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
3514             IF ( TRIM( var ) == 'ghf_eb*')  unit = 'W/m2'
3515             IF ( TRIM( var ) == 'lai*'   )  unit = 'none'
3516             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/kg*m'
3517             IF ( TRIM( var ) == 'pra*'   )  unit = 'mm'
3518             IF ( TRIM( var ) == 'prr*'   )  unit = 'mm/s'
3519             IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
3520             IF ( TRIM( var ) == 'qsws_eb*'      ) unit = 'W/m2'
3521             IF ( TRIM( var ) == 'qsws_liq_eb*'  ) unit = 'W/m2'
3522             IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
3523             IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
3524             IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'   
3525             IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
3526             IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = '' 
3527             IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = '' 
3528             IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = '' 
3529             IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
3530             IF ( TRIM( var ) == 'r_s*')     unit = 's/m' 
3531             IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
3532             IF ( TRIM( var ) == 'shf_eb*')  unit = 'W/m2'
3533             IF ( TRIM( var ) == 't*'     )  unit = 'K'
3534             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
3535             IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
3536             IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
3537
3538
3539          CASE ( 'p', 'pt', 'u', 'v', 'w' )
3540             IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
3541             IF ( TRIM( var ) == 'pt' )  unit = 'K'
3542             IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
3543             IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
3544             IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
3545             CONTINUE
3546
3547          CASE DEFAULT
3548
3549             CALL user_check_data_output( var, unit )
3550
3551             IF ( unit == 'illegal' )  THEN
3552                IF ( data_output_user(1) /= ' ' )  THEN
3553                   message_string = 'illegal value for data_output or ' //     &
3554                         'data_output_user = "' // TRIM( data_output(i) ) // '"'
3555                   CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
3556                ELSE
3557                   message_string = 'illegal value for data_output = "' //     &
3558                                    TRIM( data_output(i) ) // '"'
3559                   CALL message( 'check_parameters', 'PA0115', 1, 2, 0, 6, 0 )
3560                ENDIF
3561             ENDIF
3562
3563       END SELECT
3564!
3565!--    Set the internal steering parameters appropriately
3566       IF ( k == 0 )  THEN
3567          do3d_no(j)              = do3d_no(j) + 1
3568          do3d(j,do3d_no(j))      = data_output(i)
3569          do3d_unit(j,do3d_no(j)) = unit
3570       ELSE
3571          do2d_no(j)              = do2d_no(j) + 1
3572          do2d(j,do2d_no(j))      = data_output(i)
3573          do2d_unit(j,do2d_no(j)) = unit
3574          IF ( data_output(i)(ilen-2:ilen) == '_xy' )  THEN
3575             data_output_xy(j) = .TRUE.
3576          ENDIF
3577          IF ( data_output(i)(ilen-2:ilen) == '_xz' )  THEN
3578             data_output_xz(j) = .TRUE.
3579          ENDIF
3580          IF ( data_output(i)(ilen-2:ilen) == '_yz' )  THEN
3581             data_output_yz(j) = .TRUE.
3582          ENDIF
3583       ENDIF
3584
3585       IF ( j == 1 )  THEN
3586!
3587!--       Check, if variable is already subject to averaging
3588          found = .FALSE.
3589          DO  k = 1, doav_n
3590             IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
3591          ENDDO
3592
3593          IF ( .NOT. found )  THEN
3594             doav_n = doav_n + 1
3595             doav(doav_n) = var
3596          ENDIF
3597       ENDIF
3598
3599       i = i + 1
3600    ENDDO
3601
3602!
3603!-- Averaged 2d or 3d output requires that an averaging interval has been set
3604    IF ( doav_n > 0  .AND.  averaging_interval == 0.0_wp )  THEN
3605       WRITE( message_string, * )  'output of averaged quantity "',            &
3606                                   TRIM( doav(1) ), '_av" requires to set a ', &
3607                                   'non-zero & averaging interval'
3608       CALL message( 'check_parameters', 'PA0323', 1, 2, 0, 6, 0 )
3609    ENDIF
3610
3611!
3612!-- Check sectional planes and store them in one shared array
3613    IF ( ANY( section_xy > nz + 1 ) )  THEN
3614       WRITE( message_string, * )  'section_xy must be <= nz + 1 = ', nz + 1
3615       CALL message( 'check_parameters', 'PA0319', 1, 2, 0, 6, 0 )
3616    ENDIF
3617    IF ( ANY( section_xz > ny + 1 ) )  THEN
3618       WRITE( message_string, * )  'section_xz must be <= ny + 1 = ', ny + 1
3619       CALL message( 'check_parameters', 'PA0320', 1, 2, 0, 6, 0 )
3620    ENDIF
3621    IF ( ANY( section_yz > nx + 1 ) )  THEN
3622       WRITE( message_string, * )  'section_yz must be <= nx + 1 = ', nx + 1
3623       CALL message( 'check_parameters', 'PA0321', 1, 2, 0, 6, 0 )
3624    ENDIF
3625    section(:,1) = section_xy
3626    section(:,2) = section_xz
3627    section(:,3) = section_yz
3628
3629!
3630!-- Upper plot limit for 2D vertical sections
3631    IF ( z_max_do2d == -1.0_wp )  z_max_do2d = zu(nzt)
3632    IF ( z_max_do2d < zu(nzb+1)  .OR.  z_max_do2d > zu(nzt) )  THEN
3633       WRITE( message_string, * )  'z_max_do2d = ', z_max_do2d,                &
3634                    ' must be >= ', zu(nzb+1), '(zu(nzb+1)) and <= ', zu(nzt), &
3635                    ' (zu(nzt))'
3636       CALL message( 'check_parameters', 'PA0116', 1, 2, 0, 6, 0 )
3637    ENDIF
3638
3639!
3640!-- Upper plot limit for 3D arrays
3641    IF ( nz_do3d == -9999 )  nz_do3d = nzt + 1
3642
3643!
3644!-- Set output format string (used in header)
3645    SELECT CASE ( netcdf_data_format )
3646       CASE ( 1 )
3647          output_format_netcdf = 'netCDF classic'
3648       CASE ( 2 )
3649          output_format_netcdf = 'netCDF 64bit offset'
3650       CASE ( 3 )
3651          output_format_netcdf = 'netCDF4/HDF5'
3652       CASE ( 4 )
3653          output_format_netcdf = 'netCDF4/HDF5 classic'
3654       CASE ( 5 )
3655          output_format_netcdf = 'parallel netCDF4/HDF5'
3656       CASE ( 6 )
3657          output_format_netcdf = 'parallel netCDF4/HDF5 classic'
3658
3659    END SELECT
3660
3661#if defined( __spectra )
3662!
3663!-- Check the number of spectra level to be output
3664    i = 1
3665    DO WHILE ( comp_spectra_level(i) /= 999999  .AND.  i <= 100 )
3666       i = i + 1
3667    ENDDO
3668    i = i - 1
3669    IF ( i == 0 )  THEN
3670       WRITE( message_string, * )  'no spectra levels given'
3671       CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
3672    ENDIF
3673#endif
3674
3675!
3676!-- Check mask conditions
3677    DO mid = 1, max_masks
3678       IF ( data_output_masks(mid,1) /= ' ' .OR.                               &
3679            data_output_masks_user(mid,1) /= ' ' ) THEN
3680          masks = masks + 1
3681       ENDIF
3682    ENDDO
3683   
3684    IF ( masks < 0 .OR. masks > max_masks )  THEN
3685       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ',   &
3686            '<= ', max_masks, ' (=max_masks)'
3687       CALL message( 'check_parameters', 'PA0325', 1, 2, 0, 6, 0 )
3688    ENDIF
3689    IF ( masks > 0 )  THEN
3690       mask_scale(1) = mask_scale_x
3691       mask_scale(2) = mask_scale_y
3692       mask_scale(3) = mask_scale_z
3693       IF ( ANY( mask_scale <= 0.0_wp ) )  THEN
3694          WRITE( message_string, * )                                           &
3695               'illegal value: mask_scale_x, mask_scale_y and mask_scale_z',   &
3696               'must be > 0.0'
3697          CALL message( 'check_parameters', 'PA0326', 1, 2, 0, 6, 0 )
3698       ENDIF
3699!
3700!--    Generate masks for masked data output
3701!--    Parallel netcdf output is not tested so far for masked data, hence
3702!--    netcdf_data_format is switched back to non-paralell output.
3703       netcdf_data_format_save = netcdf_data_format
3704       IF ( netcdf_data_format > 4 )  THEN
3705          IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
3706          IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
3707          message_string = 'netCDF file formats '//                            &
3708                           '5 (parallel netCDF 4) and ' //                     &
3709                           '6 (parallel netCDF 4 Classic model) '//            &
3710                           '&are currently not supported (not yet tested) ' // &
3711                           'for masked data.&Using respective non-parallel' // & 
3712                           ' output for masked data.'
3713          CALL message( 'check_parameters', 'PA0383', 0, 0, 0, 6, 0 )
3714       ENDIF
3715       CALL init_masks
3716       netcdf_data_format = netcdf_data_format_save
3717    ENDIF
3718
3719!
3720!-- Check the NetCDF data format
3721#if ! defined ( __check )
3722    IF ( netcdf_data_format > 2 )  THEN
3723#if defined( __netcdf4 )
3724       CONTINUE
3725#else
3726       message_string = 'netCDF: netCDF4 format requested but no ' //          &
3727                        'cpp-directive __netcdf4 given & switch '  //          &
3728                        'back to 64-bit offset format'
3729       CALL message( 'check_parameters', 'PA0171', 0, 1, 0, 6, 0 )
3730       netcdf_data_format = 2
3731#endif
3732    ENDIF
3733    IF ( netcdf_data_format > 4 )  THEN
3734#if defined( __netcdf4 ) && defined( __netcdf4_parallel )
3735       CONTINUE
3736#else
3737       message_string = 'netCDF: netCDF4 parallel output requested but no ' // &
3738                        'cpp-directive __netcdf4_parallel given & switch '  // &
3739                        'back to netCDF4 non-parallel output'
3740       CALL message( 'check_parameters', 'PA0099', 0, 1, 0, 6, 0 )
3741       netcdf_data_format = netcdf_data_format - 2
3742#endif
3743    ENDIF
3744#endif
3745
3746!
3747!-- Calculate fixed number of output time levels for parallel netcdf output.
3748!-- The time dimension has to be defined as limited for parallel output,
3749!-- because otherwise the I/O performance drops significantly.
3750    IF ( netcdf_data_format > 4 )  THEN
3751
3752       ntdim_3d(0)    = INT( ( end_time - skip_time_do3d ) / dt_do3d )
3753       IF ( do3d_at_begin ) ntdim_3d(0) = ntdim_3d(0) + 1
3754       ntdim_3d(1)    = INT( ( end_time - skip_time_data_output_av )           &
3755                             / dt_data_output_av )
3756       ntdim_2d_xy(0) = INT( ( end_time - skip_time_do2d_xy ) / dt_do2d_xy )
3757       ntdim_2d_xz(0) = INT( ( end_time - skip_time_do2d_xz ) / dt_do2d_xz )
3758       ntdim_2d_yz(0) = INT( ( end_time - skip_time_do2d_yz ) / dt_do2d_yz )
3759       IF ( do2d_at_begin )  THEN
3760          ntdim_2d_xy(0) = ntdim_2d_xy(0) + 1
3761          ntdim_2d_xz(0) = ntdim_2d_xz(0) + 1
3762          ntdim_2d_yz(0) = ntdim_2d_yz(0) + 1
3763       ENDIF
3764       ntdim_2d_xy(1) = ntdim_3d(1)
3765       ntdim_2d_xz(1) = ntdim_3d(1)
3766       ntdim_2d_yz(1) = ntdim_3d(1)
3767
3768    ENDIF
3769
3770#if ! defined( __check )
3771!
3772!-- Check netcdf precison
3773    ldum = .FALSE.
3774    CALL define_netcdf_header( 'ch', ldum, 0 )
3775#endif
3776!
3777!-- Check, whether a constant diffusion coefficient shall be used
3778    IF ( km_constant /= -1.0_wp )  THEN
3779       IF ( km_constant < 0.0_wp )  THEN
3780          WRITE( message_string, * )  'km_constant = ', km_constant, ' < 0.0'
3781          CALL message( 'check_parameters', 'PA0121', 1, 2, 0, 6, 0 )
3782       ELSE
3783          IF ( prandtl_number < 0.0_wp )  THEN
3784             WRITE( message_string, * )  'prandtl_number = ', prandtl_number,  &
3785                                         ' < 0.0'
3786             CALL message( 'check_parameters', 'PA0122', 1, 2, 0, 6, 0 )
3787          ENDIF
3788          constant_diffusion = .TRUE.
3789
3790          IF ( prandtl_layer )  THEN
3791             message_string = 'prandtl_layer is not allowed with fixed ' //    &
3792                              'value of km'
3793             CALL message( 'check_parameters', 'PA0123', 1, 2, 0, 6, 0 )
3794          ENDIF
3795       ENDIF
3796    ENDIF
3797
3798!
3799!-- In case of non-cyclic lateral boundaries and a damping layer for the
3800!-- potential temperature, check the width of the damping layer
3801    IF ( bc_lr /= 'cyclic' ) THEN
3802       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3803            pt_damping_width > REAL( nx * dx ) )  THEN
3804          message_string = 'pt_damping_width out of range'
3805          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3806       ENDIF
3807    ENDIF
3808
3809    IF ( bc_ns /= 'cyclic' )  THEN
3810       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3811            pt_damping_width > REAL( ny * dy ) )  THEN
3812          message_string = 'pt_damping_width out of range'
3813          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3814       ENDIF
3815    ENDIF
3816
3817!
3818!-- Check value range for rif
3819    IF ( rif_min >= rif_max )  THEN
3820       WRITE( message_string, * )  'rif_min = ', rif_min, ' must be less ',    &
3821                                   'than rif_max = ', rif_max
3822       CALL message( 'check_parameters', 'PA0125', 1, 2, 0, 6, 0 )
3823    ENDIF
3824
3825!
3826!-- Check random generator
3827    IF ( (random_generator /= 'system-specific'     .AND.                      &
3828          random_generator /= 'random-parallel'   ) .AND.                      &
3829          random_generator /= 'numerical-recipes' )  THEN
3830       message_string = 'unknown random generator: random_generator = "' //    &
3831                        TRIM( random_generator ) // '"'
3832       CALL message( 'check_parameters', 'PA0135', 1, 2, 0, 6, 0 )
3833    ENDIF
3834
3835!
3836!-- Determine upper and lower hight level indices for random perturbations
3837    IF ( disturbance_level_b == -9999999.9_wp )  THEN
3838       IF ( ocean ) THEN
3839          disturbance_level_b     = zu((nzt*2)/3)
3840          disturbance_level_ind_b = ( nzt * 2 ) / 3
3841       ELSE
3842          disturbance_level_b     = zu(nzb+3)
3843          disturbance_level_ind_b = nzb + 3
3844       ENDIF
3845    ELSEIF ( disturbance_level_b < zu(3) )  THEN
3846       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3847                           disturbance_level_b, ' must be >= ', zu(3), '(zu(3))'
3848       CALL message( 'check_parameters', 'PA0126', 1, 2, 0, 6, 0 )
3849    ELSEIF ( disturbance_level_b > zu(nzt-2) )  THEN
3850       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3851                   disturbance_level_b, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3852       CALL message( 'check_parameters', 'PA0127', 1, 2, 0, 6, 0 )
3853    ELSE
3854       DO  k = 3, nzt-2
3855          IF ( disturbance_level_b <= zu(k) )  THEN
3856             disturbance_level_ind_b = k
3857             EXIT
3858          ENDIF
3859       ENDDO
3860    ENDIF
3861
3862    IF ( disturbance_level_t == -9999999.9_wp )  THEN
3863       IF ( ocean )  THEN
3864          disturbance_level_t     = zu(nzt-3)
3865          disturbance_level_ind_t = nzt - 3
3866       ELSE
3867          disturbance_level_t     = zu(nzt/3)
3868          disturbance_level_ind_t = nzt / 3
3869       ENDIF
3870    ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
3871       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3872                   disturbance_level_t, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3873       CALL message( 'check_parameters', 'PA0128', 1, 2, 0, 6, 0 )
3874    ELSEIF ( disturbance_level_t < disturbance_level_b )  THEN
3875       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3876                   disturbance_level_t, ' must be >= disturbance_level_b = ',  &
3877                   disturbance_level_b
3878       CALL message( 'check_parameters', 'PA0129', 1, 2, 0, 6, 0 )
3879    ELSE
3880       DO  k = 3, nzt-2
3881          IF ( disturbance_level_t <= zu(k) )  THEN
3882             disturbance_level_ind_t = k
3883             EXIT
3884          ENDIF
3885       ENDDO
3886    ENDIF
3887
3888!
3889!-- Check again whether the levels determined this way are ok.
3890!-- Error may occur at automatic determination and too few grid points in
3891!-- z-direction.
3892    IF ( disturbance_level_ind_t < disturbance_level_ind_b )  THEN
3893       WRITE( message_string, * )  'disturbance_level_ind_t = ',               &
3894                disturbance_level_ind_t, ' must be >= disturbance_level_b = ', &
3895                disturbance_level_b
3896       CALL message( 'check_parameters', 'PA0130', 1, 2, 0, 6, 0 )
3897    ENDIF
3898
3899!
3900!-- Determine the horizontal index range for random perturbations.
3901!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed
3902!-- near the inflow and the perturbation area is further limited to ...(1)
3903!-- after the initial phase of the flow.
3904   
3905    IF ( bc_lr /= 'cyclic' )  THEN
3906       IF ( inflow_disturbance_begin == -1 )  THEN
3907          inflow_disturbance_begin = MIN( 10, nx/2 )
3908       ENDIF
3909       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )&
3910       THEN
3911          message_string = 'inflow_disturbance_begin out of range'
3912          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3913       ENDIF
3914       IF ( inflow_disturbance_end == -1 )  THEN
3915          inflow_disturbance_end = MIN( 100, 3*nx/4 )
3916       ENDIF
3917       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )    &
3918       THEN
3919          message_string = 'inflow_disturbance_end out of range'
3920          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3921       ENDIF
3922    ELSEIF ( bc_ns /= 'cyclic' )  THEN
3923       IF ( inflow_disturbance_begin == -1 )  THEN
3924          inflow_disturbance_begin = MIN( 10, ny/2 )
3925       ENDIF
3926       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )&
3927       THEN
3928          message_string = 'inflow_disturbance_begin out of range'
3929          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3930       ENDIF
3931       IF ( inflow_disturbance_end == -1 )  THEN
3932          inflow_disturbance_end = MIN( 100, 3*ny/4 )
3933       ENDIF
3934       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )    &
3935       THEN
3936          message_string = 'inflow_disturbance_end out of range'
3937          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3938       ENDIF
3939    ENDIF
3940
3941    IF ( random_generator == 'random-parallel' )  THEN
3942       dist_nxl = nxl;  dist_nxr = nxr
3943       dist_nys = nys;  dist_nyn = nyn
3944       IF ( bc_lr == 'radiation/dirichlet' )  THEN
3945          dist_nxr    = MIN( nx - inflow_disturbance_begin, nxr )
3946          dist_nxl(1) = MAX( nx - inflow_disturbance_end, nxl )
3947       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
3948          dist_nxl    = MAX( inflow_disturbance_begin, nxl )
3949          dist_nxr(1) = MIN( inflow_disturbance_end, nxr )
3950       ENDIF
3951       IF ( bc_ns == 'dirichlet/radiation' )  THEN
3952          dist_nyn    = MIN( ny - inflow_disturbance_begin, nyn )
3953          dist_nys(1) = MAX( ny - inflow_disturbance_end, nys )
3954       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
3955          dist_nys    = MAX( inflow_disturbance_begin, nys )
3956          dist_nyn(1) = MIN( inflow_disturbance_end, nyn )
3957       ENDIF
3958    ELSE
3959       dist_nxl = 0;  dist_nxr = nx
3960       dist_nys = 0;  dist_nyn = ny
3961       IF ( bc_lr == 'radiation/dirichlet' )  THEN
3962          dist_nxr    = nx - inflow_disturbance_begin
3963          dist_nxl(1) = nx - inflow_disturbance_end
3964       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
3965          dist_nxl    = inflow_disturbance_begin
3966          dist_nxr(1) = inflow_disturbance_end
3967       ENDIF
3968       IF ( bc_ns == 'dirichlet/radiation' )  THEN
3969          dist_nyn    = ny - inflow_disturbance_begin
3970          dist_nys(1) = ny - inflow_disturbance_end
3971       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
3972          dist_nys    = inflow_disturbance_begin
3973          dist_nyn(1) = inflow_disturbance_end
3974       ENDIF
3975    ENDIF
3976
3977!
3978!-- A turbulent inflow requires Dirichlet conditions at the respective inflow
3979!-- boundary (so far, a turbulent inflow is realized from the left side only)
3980    IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' )  THEN
3981       message_string = 'turbulent_inflow = .T. requires a Dirichlet ' //      &
3982                        'condition at the inflow boundary'
3983       CALL message( 'check_parameters', 'PA0133', 1, 2, 0, 6, 0 )
3984    ENDIF
3985
3986!
3987!-- Turbulent inflow requires that 3d arrays have been cyclically filled with
3988!-- data from prerun in the first main run
3989    IF ( turbulent_inflow  .AND.  initializing_actions /= 'cyclic_fill'  .AND. &
3990         initializing_actions /= 'read_restart_data' )  THEN
3991       message_string = 'turbulent_inflow = .T. requires ' //                  &
3992                        'initializing_actions = ''cyclic_fill'' '
3993       CALL message( 'check_parameters', 'PA0055', 1, 2, 0, 6, 0 )
3994    ENDIF
3995
3996!
3997!-- In case of turbulent inflow calculate the index of the recycling plane
3998    IF ( turbulent_inflow )  THEN
3999       IF ( recycling_width == 9999999.9_wp )  THEN
4000!
4001!--       Set the default value for the width of the recycling domain
4002          recycling_width = 0.1_wp * nx * dx
4003       ELSE
4004          IF ( recycling_width < dx  .OR.  recycling_width > nx * dx )  THEN
4005             WRITE( message_string, * )  'illegal value for recycling_width:', &
4006                                         ' ', recycling_width
4007             CALL message( 'check_parameters', 'PA0134', 1, 2, 0, 6, 0 )
4008          ENDIF
4009       ENDIF
4010!
4011!--    Calculate the index
4012       recycling_plane = recycling_width / dx
4013    ENDIF
4014
4015!
4016!-- Determine damping level index for 1D model
4017    IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
4018       IF ( damp_level_1d == -1.0_wp )  THEN
4019          damp_level_1d     = zu(nzt+1)
4020          damp_level_ind_1d = nzt + 1
4021       ELSEIF ( damp_level_1d < 0.0_wp  .OR.  damp_level_1d > zu(nzt+1) )  THEN
4022          WRITE( message_string, * )  'damp_level_1d = ', damp_level_1d,       &
4023                 ' must be > 0.0 and < ', zu(nzt+1), '(zu(nzt+1))'
4024          CALL message( 'check_parameters', 'PA0136', 1, 2, 0, 6, 0 )
4025       ELSE
4026          DO  k = 1, nzt+1
4027             IF ( damp_level_1d <= zu(k) )  THEN
4028                damp_level_ind_1d = k
4029                EXIT
4030             ENDIF
4031          ENDDO
4032       ENDIF
4033    ENDIF
4034
4035!
4036!-- Check some other 1d-model parameters
4037    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                   &
4038         TRIM( mixing_length_1d ) /= 'blackadar' )  THEN
4039       message_string = 'mixing_length_1d = "' // TRIM( mixing_length_1d ) //  &
4040                        '" is unknown'
4041       CALL message( 'check_parameters', 'PA0137', 1, 2, 0, 6, 0 )
4042    ENDIF
4043    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND.                     &
4044         TRIM( dissipation_1d ) /= 'detering' )  THEN
4045       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //      &
4046                        '" is unknown'
4047       CALL message( 'check_parameters', 'PA0138', 1, 2, 0, 6, 0 )
4048    ENDIF
4049
4050!
4051!-- Set time for the next user defined restart (time_restart is the
4052!-- internal parameter for steering restart events)
4053    IF ( restart_time /= 9999999.9_wp )  THEN
4054       IF ( restart_time > time_since_reference_point )  THEN
4055          time_restart = restart_time
4056       ENDIF
4057    ELSE
4058!
4059!--    In case of a restart run, set internal parameter to default (no restart)
4060!--    if the NAMELIST-parameter restart_time is omitted
4061       time_restart = 9999999.9_wp
4062    ENDIF
4063
4064!
4065!-- Set default value of the time needed to terminate a model run
4066    IF ( termination_time_needed == -1.0_wp )  THEN
4067       IF ( host(1:3) == 'ibm' )  THEN
4068          termination_time_needed = 300.0_wp
4069       ELSE
4070          termination_time_needed = 35.0_wp
4071       ENDIF
4072    ENDIF
4073
4074!
4075!-- Check the time needed to terminate a model run
4076    IF ( host(1:3) == 't3e' )  THEN
4077!
4078!--    Time needed must be at least 30 seconds on all CRAY machines, because
4079!--    MPP_TREMAIN gives the remaining CPU time only in steps of 30 seconds
4080       IF ( termination_time_needed <= 30.0_wp )  THEN
4081          WRITE( message_string, * )  'termination_time_needed = ',            &
4082                 termination_time_needed, ' must be > 30.0 on host "',         &
4083                 TRIM( host ), '"'
4084          CALL message( 'check_parameters', 'PA0139', 1, 2, 0, 6, 0 )
4085       ENDIF
4086    ELSEIF ( host(1:3) == 'ibm' )  THEN
4087!
4088!--    On IBM-regatta machines the time should be at least 300 seconds,
4089!--    because the job time consumed before executing palm (for compiling,
4090!--    copying of files, etc.) has to be regarded
4091       IF ( termination_time_needed < 300.0_wp )  THEN
4092          WRITE( message_string, * )  'termination_time_needed = ',            &
4093                 termination_time_needed, ' should be >= 300.0 on host "',     &
4094                 TRIM( host ), '"'
4095          CALL message( 'check_parameters', 'PA0140', 1, 2, 0, 6, 0 )
4096       ENDIF
4097    ENDIF
4098
4099!
4100!-- Check pressure gradient conditions
4101    IF ( dp_external .AND. conserve_volume_flow )  THEN
4102       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo', &
4103            'w are .TRUE. but one of them must be .FALSE.'
4104       CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
4105    ENDIF
4106    IF ( dp_external )  THEN
4107       IF ( dp_level_b < zu(nzb) .OR. dp_level_b > zu(nzt) )  THEN
4108          WRITE( message_string, * )  'dp_level_b = ', dp_level_b, ' is out ', &
4109               ' of range'
4110          CALL message( 'check_parameters', 'PA0151', 1, 2, 0, 6, 0 )
4111       ENDIF
4112       IF ( .NOT. ANY( dpdxy /= 0.0_wp ) )  THEN
4113          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze', &
4114               'ro, i.e. the external pressure gradient & will not be applied'
4115          CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
4116       ENDIF
4117    ENDIF
4118    IF ( ANY( dpdxy /= 0.0_wp ) .AND. .NOT. dp_external )  THEN
4119       WRITE( message_string, * )  'dpdxy is nonzero but dp_external is ',     &
4120            '.FALSE., i.e. the external pressure gradient & will not be applied'
4121       CALL message( 'check_parameters', 'PA0153', 0, 1, 0, 6, 0 )
4122    ENDIF
4123    IF ( conserve_volume_flow )  THEN
4124       IF ( TRIM( conserve_volume_flow_mode ) == 'default' )  THEN
4125
4126          conserve_volume_flow_mode = 'initial_profiles'
4127
4128       ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles' .AND.  &
4129            TRIM( conserve_volume_flow_mode ) /= 'inflow_profile' .AND.        &
4130            TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' )  THEN
4131          WRITE( message_string, * )  'unknown conserve_volume_flow_mode: ',   &
4132               conserve_volume_flow_mode
4133          CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 )
4134       ENDIF
4135       IF ( (bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic')  .AND.                &
4136          TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
4137          WRITE( message_string, * )  'non-cyclic boundary conditions ',       &
4138               'require  conserve_volume_flow_mode = ''initial_profiles'''
4139          CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 )
4140       ENDIF
4141       IF ( bc_lr == 'cyclic'  .AND.  bc_ns == 'cyclic'  .AND.                 &
4142            TRIM( conserve_volume_flow_mode ) == 'inflow_profile' )  THEN
4143          WRITE( message_string, * )  'cyclic boundary conditions ',           &
4144               'require conserve_volume_flow_mode = ''initial_profiles''',     &
4145               ' or ''bulk_velocity'''
4146          CALL message( 'check_parameters', 'PA0156', 1, 2, 0, 6, 0 )
4147       ENDIF
4148    ENDIF
4149    IF ( ( u_bulk /= 0.0_wp .OR. v_bulk /= 0.0_wp ) .AND.                      &
4150         ( .NOT. conserve_volume_flow .OR.                                     &
4151         TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )  THEN
4152       WRITE( message_string, * )  'nonzero bulk velocity requires ',          &
4153            'conserve_volume_flow = .T. and ',                                 &
4154            'conserve_volume_flow_mode = ''bulk_velocity'''
4155       CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
4156    ENDIF
4157
4158!
4159!-- Check particle attributes
4160    IF ( particle_color /= 'none' )  THEN
4161       IF ( particle_color /= 'absuv'  .AND.  particle_color /= 'pt*'  .AND.   &
4162            particle_color /= 'z' )  THEN
4163          message_string = 'illegal value for parameter particle_color: ' //   &
4164                           TRIM( particle_color)
4165          CALL message( 'check_parameters', 'PA0313', 1, 2, 0, 6, 0 )
4166       ELSE
4167          IF ( color_interval(2) <= color_interval(1) )  THEN
4168             message_string = 'color_interval(2) <= color_interval(1)'
4169             CALL message( 'check_parameters', 'PA0315', 1, 2, 0, 6, 0 )
4170          ENDIF
4171       ENDIF
4172    ENDIF
4173
4174    IF ( particle_dvrpsize /= 'none' )  THEN
4175       IF ( particle_dvrpsize /= 'absw' )  THEN
4176          message_string = 'illegal value for parameter particle_dvrpsize:' // &
4177                           ' ' // TRIM( particle_color)
4178          CALL message( 'check_parameters', 'PA0314', 1, 2, 0, 6, 0 )
4179       ELSE
4180          IF ( dvrpsize_interval(2) <= dvrpsize_interval(1) )  THEN
4181             message_string = 'dvrpsize_interval(2) <= dvrpsize_interval(1)'
4182             CALL message( 'check_parameters', 'PA0316', 1, 2, 0, 6, 0 )
4183          ENDIF
4184       ENDIF
4185    ENDIF
4186
4187!
4188!-- Check nudging and large scale forcing from external file
4189    IF ( nudging .AND. ( .NOT. large_scale_forcing ) )  THEN
4190       message_string = 'Nudging requires large_scale_forcing = .T.. &'//      &
4191                        'Surface fluxes and geostrophic wind should be &'//    &
4192                        'prescribed in file LSF_DATA'
4193       CALL message( 'check_parameters', 'PA0374', 1, 2, 0, 6, 0 )
4194    ENDIF
4195
4196    IF ( large_scale_forcing .AND. ( bc_lr /= 'cyclic'  .OR.                   &
4197                                    bc_ns /= 'cyclic' ) )  THEN
4198       message_string = 'Non-cyclic lateral boundaries do not allow for &' //  &
4199                        'the usage of large scale forcing from external file.'
4200       CALL message( 'check_parameters', 'PA0375', 1, 2