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

Last change on this file since 3294 was 3294, checked in by raasch, 6 years ago

modularization of the ocean code

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