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

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

Bugfix in large-scale forcing

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