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

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

various changes to avoid compiler warnings (mainly removal of unused variables)

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