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

Last change on this file since 2332 was 2329, checked in by knoop, 7 years ago

Bugfix for topography usage with anelastic approximation and boussinesq approximation with air density != 1

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