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

Last change on this file since 2897 was 2883, checked in by Giersch, 7 years ago

Format 325 of header output has been changed, dt_dopr_listing is not set to the default value zero anymore

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