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

Last change on this file since 3301 was 3301, checked in by gronemeier, 6 years ago

corrected former revision section

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