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

Last change on this file since 2675 was 2669, checked in by raasch, 7 years ago

file attributes and activation strings in .palm.iofiles revised, file extensions for nesting, masked output, wind turbine data, etc. changed

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