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

Last change on this file since 3582 was 3582, checked in by suehring, 5 years ago

Merge branch salsa with trunk

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