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

Last change on this file since 3619 was 3597, checked in by maronga, 6 years ago

revised calculation of near surface air potential temperature

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