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

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

Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled

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