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

Last change on this file since 1765 was 1765, checked in by raasch, 7 years ago

last commit documented

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