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

Last change on this file since 2918 was 2918, checked in by gronemeier, 3 years ago

moved initialization of mixing length to turbulence_closure_mod; consider walls in calculation of ml in rans-mode

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