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

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

merge with branch rans: update of rans mode and data output

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