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

Last change on this file since 2550 was 2550, checked in by boeske, 7 years ago

enable simulations with complex terrain

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