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

Last change on this file since 2575 was 2575, checked in by maronga, 6 years ago

bugfix in radiation model and improvements in land surface scheme

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