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

Last change on this file since 2934 was 2934, checked in by suehring, 3 years ago

Synchronize parent and child model after initialization and spinup phase; Check for consistent setting of spinup times in parent and child model; remove obsolete masking of tendency arrays during initialization

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