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

Last change on this file since 2743 was 2743, checked in by suehring, 6 years ago

In case of natural- and urban-type surfaces output surfaces fluxes in W/m2

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