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

Last change on this file since 2504 was 2422, checked in by raasch, 7 years ago

new palm-scripts set to version number 1.0, further check for consistent restart settings

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