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

Last change on this file since 2513 was 2513, checked in by kanani, 7 years ago

Bugfixes related to profile output of passive scalar

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