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

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

Error message PA0476 added. Bugfix in case of coupled runs.

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