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

Last change on this file since 2354 was 2354, checked in by schwenkel, 7 years ago

bugfix for lsm data output

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