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

Last change on this file since 2339 was 2339, checked in by gronemeier, 7 years ago

corrected timestamp in header

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