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

Last change on this file since 3034 was 3034, checked in by raasch, 6 years ago

bugfix: check that initializing_actions has been set

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