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

Last change on this file since 2512 was 2508, checked in by suehring, 7 years ago

Bugfixes in SGS-TKE buoyancy production; revised initialization of vertical-gradient levels in case of ocean runs

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