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

Last change on this file since 2410 was 2375, checked in by schwenkel, 7 years ago

improved aerosol initialization for bulk microphysics

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