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

Last change on this file since 3039 was 3035, checked in by schwenkel, 7 years ago

Add option to initialize warm air bubble close to surface

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