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

Last change on this file since 3525 was 3525, checked in by kanani, 3 years ago

Changes related to clean-up of biometeorology (by dom_dwd_user)

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