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

Last change on this file since 2938 was 2938, checked in by suehring, 4 years ago

Nesting in RANS-LES and RANS-RANS mode enabled; synthetic turbulence generator at all lateral boundaries in nesting or non-cyclic forcing mode; revised Inifor initialization in nesting mode

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