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

Last change on this file since 3082 was 3065, checked in by Giersch, 6 years ago

New vertical stretching procedure has been introduced

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