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

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

add variable description

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