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

Last change on this file since 2085 was 2085, checked in by knoop, 5 years ago

last commit documented

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