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

Last change on this file since 2792 was 2776, checked in by Giersch, 6 years ago

Skipping of module related restart data changed + adapting synthetic turbulence generator to current restart procedure

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