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

Last change on this file since 1691 was 1691, checked in by maronga, 9 years ago

various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings

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