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

Last change on this file since 3575 was 3575, checked in by kanani, 5 years ago

Bugfix for output time levels in parallel NetCDF output with spinup (check_parameters)

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