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

Last change on this file since 3232 was 3232, checked in by raasch, 6 years ago

references to mrun replaced by palmrun, and updated

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