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

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

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

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