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

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

last commit documented

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