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

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

Fix for biomet output (ticket:757), merge of uv_exposure into biometeorology_mod

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