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

Last change on this file since 3026 was 2980, checked in by suehring, 6 years ago

Improvement for spinup checks concering consistency in nesting mode

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