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

Last change on this file since 2663 was 2628, checked in by schwenkel, 6 years ago

enable particle advection with grid stretching and some formatation changes in lpm

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