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

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

Remaining error messages revised, comments extended

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