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

Last change on this file since 2821 was 2817, checked in by knoop, 6 years ago

Preliminary gust module interface implemented

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