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

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

Bugfix in if query.

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