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

Last change on this file since 2320 was 2320, checked in by suehring, 4 years ago

large-scale forcing and nudging modularized

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