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

Last change on this file since 3593 was 3589, checked in by suehring, 6 years ago

Remove erroneous UTF encoding; last commit documented

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